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Models of galaxy and halo clustering commonly assume that the tracers can be treated as a 
continuous field locally biased with respect to the underlying mass distribution. In the peak model 
pioneered by Bardeen et al. (1986), one considers instead density maxima of the initial, Gaussian 
mass density field as an approximation to the formation site of virialized objects. In this paper, 
CNj ■ the peak model is extended in two ways to improve its predictive accuracy. Firstly, we derive the 

two-point correlation function of initial density peaks up to second order and demonstrate that a 
O , peak-background split approach can be applied to obtain the fc-independent and fc-dependent peak 

bias factors at all orders. Secondly, we explore the gravitational evolution of the peak correlation 
function within the Zel'dovich approximation. We show that the local (Lagrangian) bias approach 
emerges as a special case of the peak model, in which all bias parameters are scale-independent and 
there is no statistical velocity bias. We apply our formulae to study how the Lagrangian peak biasing, 
the diffusion due to large scale flows and the mode-coupling due to nonlocal interactions affect the 
scale dependence of bias from small separations up to the baryon acoustic oscillation (BAO) scale. 
For 2a density peaks collapsing at z = 0.3, our model predicts a ~5% residual scale-dependent bias 
around the acoustic scale that arises mostly from first-order Lagrangian peak biasing (as opposed 
pH ' to second-order gravity mode-coupling). We also search for a scale dependence of bias in the large 

Q scale auto-correlation of massive halos extracted from a very large N-body simulation provided by 

the MICE collaboration. For halos with mass M > 10 14 Mq//i, our measurements demonstrate a 
scale-dependent bias across the BAO feature which is very well reproduced by a prediction based 
on the peak model. 
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PACS numbers: 98.80.-k, 98.65.-r, 95.35. +d, 98.80.Es 



I. INTRODUCTION 



A considerable amount of effort has already been invested in measuring the large scale distribution of galaxies, 
0^ ' especially the galaxy two-point correlation function and power spectrum, to constrain viable cosmological models 
[e.g., [ll-[l4j. The amplitude, shape and baryon acoustic feature in these two-point statistics encode a wealth of 
cosmological information [15l432l |. Ongoing and planned galaxy surveys of the high redshift Universe will furnish 
measurements of the underlying mass distribution with unprecedent precision and statistics. Alongside this great 
observational effort, interpreting this vast amount of data will require a much better understanding of the relation 
between the surveyed galaxies and the mass fluctuations they are thought to trace. 

Essentially all models of galaxy clustering assume that galaxies are biased tracers of the mass density fluctuation 
field. Although this bias is expected to be nonlinear, scale-dependent and stochastic [H,[34j|, the si mpler, li near, scale- 
independent, deterministic model has proved to be an extremely useful first order approximation |35l437j |. However, 
in order to predict corrections beyond linear order to the galaxy two-point correlation, or even the leading order 
contribution to higher order statistics such as the galaxy three-point correlation or bispectrum, one must address the 
complications which arise from nonlinearity, scale dependence and stochasticity. For example, if the bias relation is 
established in coordinate space, then nonlinear biasing will produce scale dependence and stochasticity in Fourier 
space, and vice- versa [e.g., [38|. This randomness will add to other sources of stochasticity which may arise, for 
example, from the fact that the formation of galaxies and halos depends on quantities other than the mass density 
field (e.g., the surrounding tidal field). Moreover, the bias may be established in the initial conditions (Lagrangian 
bias) or, alternatively, at the present time (Eulerian bias). In the former case, the bias between the tracers and the 
mass will be affected by the subsequent, nonlinear gravitational evolution. This will introduce additional nonlinearity, 
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scale dependence and stochasticity. Furthermore, if the velocities of the tracers differ from those of the mass elements, 
then this will complicate the application of the continuity equation to describe the redshift evolution of bias. 

Current analytic approaches to galaxy and dark matter halo clustering take into account some of these complications. 
In most models, the fundamental quantity is the overdensity of tracers S g ^(R,x) within a sphere of radius R centered 
at position x. It is commonly assumed that (5 gj h is solely a function of the local mass overdensity S m [3^, H^j [see also 
l4lj |. whose Taylor expansion coefficients are the galaxy or halo bias parameters 6jv [4^ - |. If this bias is established 
at a different time than the epoch at which the tracers are observed, then this local bias scheme is combined with 
some (Eulerian or Lagrangian) perturbative treatment of gravitational instability [see [45l. for a review of perturbation 
theory] to predict the galaxy or halo power spectrum, bispectrum etc. [e.g.. l43l . l46l - l58j . This formalism can be extended 
to include stochasticity by formulating the biasing in terms of the conditional probability distribution P(<5 gi h|<5 m ) of 
<5 g .h at a given <5 m [e.g., [HI, H^, [59l - [61| . One of the main caveats with such local biasing schemes is that galaxies 
(or halos) are treated as though they define a continuous field smoothed on some scale R, whereas they are, in fact, 
discrete objects. 

The peaks approach to galaxy and dark matter halo clustering is interesting because it exhibits all of the compli- 
cations mentioned above while also accounting for the discrete nature of the tracers (after all, peaks define a point 
process). In this model, the fundamental quantity is the set of positions which are local maxima of the density field 
(from which a peak overabundance Sp^R, x) in spheres of radius R could in principle be derived). Since the evolved 
density field is highly nonlinear, the peak constraint is generally applied to the initial (Lagrangian) Gaussian density 
field, with the assumption that the most prominent peaks should be in one-to-one correspondence with luminous 
galaxies or massive halos in the low redshift Universe [see, e.g.. |62[ - I64l for numerical studies of this association] . Peak 
abundances, profiles and correlation functions in real and redshift space have been studied in the literature [36l. I65l475| . 
Some of these results have been used to interpret the abundance and cluste ring of rich clusters 4^, 7p-[80| , constrain 



the power spectrum of mass fluctuations [81], 182] , and study evolution bias [83j and assembly bias [84 

On asymptoti cally large scales, peaks are linearly biased tracers of the mass density field, and this bias is scale 
independent [35|, HfJ H3, ■ However, these conclusions are based on a configuration space argument known as the 
peak background split - which establishes a relation between the sensitivity of the peak bias factors and the peak 
abundances on the peak height - whereas a Fourier space analysis suggests that the linear bias factor of peaks is the 
sum of two terms, one of which is fc-dependent [38l . |74| . In configuration space, this leads to scale dependence of the 
bias and stochasticity. The fc-dependence of the linear peak bias arises from the peak constraint, i.e. the fact that one 
must specify not only the value of the mass density field but also its first two derivatives to define a peak. Therefore, 
this is a model in which the bias depends on quantities other than the local density. Moreover, as mentioned above, 
the peak biasing is applied to the initial Gaussian density field so that the late time peak bias is modified by nonlinear 
evolution and associated stochasticity. In this regards, peaks exhibit a nontrivial velocity bias [75j ]. which further 
complicates the nonlinear evolution. 

In the peak model, both the constant and the fc-dependent piece of the linear bias factor depend on peak height. 
As shown in (75[, the scale-independent contribution can be derived from the peak background split argument. In 
the first half of this paper, we demonstrate that the Fourier space approach also predicts constant and fc-dependent 
contributions to the second and higher order peak bias factors. We then show that the scale-independent parts of 
all these nonlinear peak bias factors can be also derived from the peak background split argument, thus generalizing 
the result of [75j]. We go on to show how the peak background split approach can be used to determine the scale- 
dependent part of the peak bias factors, first at linear order, and then for all nonlinear orders as well. This is 
particularly interesting because it illustrates how the peak background split argument should be implemented if the 
abundance of the biased tracers (in this case, peaks) depends on quantities other than the local mass overdensity (in 
this case, the first and second derivatives of the mass density field). 

As recognized in [74| , the fc-dependence of the first order peak bias strongly amplifies the contrast of the baryon 
acoustic oscillation [or BAO. see l86l. and references therein] in the correlation of initial density maxima. However, 
this calculation was performed for peaks identified in the initial conditions, so there was no clear connection with the 
clustering of dark matter halos and galaxies. This is also true of all the results presented in the first half of this paper. 
To remedy this problem, we show in the second half how the effects of the (nonlinear, nonlocal) gravitational evolution 
of density peaks can be incorporated in the peak model. This allows us to ascertain the extent to which the initial 
scale dependence of bias across the BAO survives at late times. Our analysis incorporates two main complications 
that are usually ignored in local bias schemes. Namely, peak biasing depends on more than just the value of the local 
density, and peaks exhibit a velocity bias which (in addition to merging) complicates analyses based on the continuity 
equation. Finally, we show that taking into account these effects is of more than academic interest: Our peaks model 
provides a very good description of the scale dependence of the bias of massive halos in numerical simulations - halos 
that are expected to host the luminous red galaxies (LRGs) which are often targeted in BAO experiments. 

The paper is organized as follows. Section ^TT] briefly reviews known results and introduce some useful definitions. 
Section i jllll focuses on the correlation of initial density peaks of a Gaussian random field. It is shown that the 
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scale-dependent and scale-independent parts of the peak bias parameters can be derived from a peak-background 
split argument. Section EjlVI considers the gravitational evolution of the peak correlation function in the Zel'dovich 
approximation. It is shown that, in addition to gravity mode-coupling, the Lagrangian peak biasing can generate 
a significant scale-dependent bias across the baryonic acoustic feature at the collapse epoch. Measurements of bias 
at BAO scales from the clustering of massive halos are also presented and compared with the model. Section fJV] 
summarizes our results. Technical details of the calculation can be found in Appendix ElAlandlBl 



II. DEFINITIONS, NOTATIONS AND KNOWN RESULTS 

We begin by introducing some definitions and reviewing known results about the clustering of density peaks in 
Gaussian random fields. Next, we derive the peak correlation at second order. This result will serve as input to the 
calculation of the evolved correlation of density peaks. 



A. Spectral moments 



The statistical properties of density peaks depend not only on the underlying density field, but also on its first 
and second derivatives. We are, therefore, interested in the linear (Gaussian) density field S and its first and second 
derivatives, diS and didjS. In this regard, it is convenient to introduce the normalized variables v = 8/<jq, r\i = diS/ai 
and Qj = didjS / where the a n are the spectral moments of the matter power spectrum, 

1 r°° 

al(R s , zq) = — g / dk fc 2 (" +1 > P s (k, z )W(kRs) 2 . (1) 
2tt 2 Jo 

P$(k,zo) denotes the dimensionless power spectrum of the linear density field at redshift zq, and W is a spherically 
symmetric smoothing kernel of length Rs introduced to ensure convergence of all spectral moments. A Gaussian filter 
will be adopted throughout this paper. We will use the notation Ps s (k, zq) to denote Ps(k, zq)W (kRs) 2 ■ The ratio 
(To/fi is proportional to the typical separation between zero-crossings of the density field [36(. For subsequent use, 
we also define the spectral parameters 

ln(Rs) = (2) 

which reflect the range over which k 2 ^ n ~^Ps s (k, zq) is large. We will also work with the scaled velocities v, = 
v P i/(aH f) and with the curvature u = —d 2 8/<72- Here, v p ,(x) is the i-th component of the (proper) peculiar velocity, 
H = dhua/dt, f = d\nD/d\na is the logarithmic derivative of the linear theory growth rate D(zq) and d 2 = d l di is the 
Laplacian. Note that v, has dimensions of length. 

The analogous quantities to er 2 at non-zero separation are defined as follows: 

$ n) (R s ,r,z ) = ^ dkk 2 ^P Ss (k,z )Mkr), (3) 

where ji{x) are spherical Bessel functions. As £ gets larger, these harmonic transforms become increasingly sensitive 
to small scale power. The auto- and cross-correlations of the fields Vi(x), ??i(x), ^(x), m(x) and Cij( x ) can generally be 
decomposed into components with definite transformation properties under rotations. [74| gives explicit expressions 
for the isotropic and homogeneous linear density field. 



B. Peak biasing and 2-point correlation function at the first order 

Although density peaks form a well-behaved point process, the large scale asympotics of the two-point correlation 
£pk( r i z) and line of sight mean streaming vi 2 {v, z) ■ f of peaks of height v and curvature u identified on a scale Rs in 
the initial Gaussian density field linearly extrapolated at redshift zq can be thought of as arising from the continuous, 
deterministic bias relation [74], [7|| 

Sp k (iy, u, R s , x) = Ms(x, z ) - fo c 9 2 <5s(x, z ), (4) 

2 

V pk (i?S, X, Zq) = V S (x, Zq) 1 cMs(x, Zq) , (5) 
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which is nonlocal owing to the smoothing of the mass distribution. Here, (5 p k and v p k are the average peak ovcrdcnsity 
and velocity, <5g and v# are the mass density and velocity smoothed at scale Rs (so as to retain only the large scale, 
coherent motion of the peak), d 2 is the Laplacian, and the bias parameters b v and b^ are 

b v (v,u,R s ,zo) = — — -(- — ) , b c (v,u,R s ,z Q ) = — — - ( - — ^) . (6) 

<t (R s , z Q ) V 1 - 7i / cr 2 (i?s, 20) V 1 - 75 tJ 

The bias coefficient b v is dimensionless, whereas b^ has units of (length) 2 . In fact, b v is precisely the amplification 
factor found by [36| who neglected derivatives of the density correlation function (i.e. their analysis assumes b^ = 0). 
Unlike v p k, <5 p k does not depend on z (as expected) because the redshift dependence of b v , b^ cancels the factor D(z ) 
coming from <5s(x, z )- Note also that, if b v > the effective peak density <5 p k(x) can be less than -1 in deep voids. 
However, this is not a problem because <5 p k(x) is not an observable quantity (this is not a count-in-cell density). 

In what follows, we will focus on the clustering of initial density peaks of significance v, for which the first order 
bias parameters are 



1 ( v — 71 u 
a {R s ,z ) V 1 - 7i 



b v {u, Rs, zo) = _ , D - ( V^TT ) (7) 



\(y, R S , z ) = —± , (^t) ■ (8) 

a 2 {Rs,z Q ) V 1 - it J 

Here, the overline denotes the averaging over the peak curvature, so that u = u{u, Rs) is the mean curvature of peaks 
of height v on filtering scale Rs- It is convenient to define the quantity 6 sp k as the Fourier space multiplication by 

&spk(*:, z ) = b v + b c k 2 , (9) 

where we have omitted the explicit dependence on is, u and Rs for brevity. Although 6 S pk(fc, zo) has the same functional 
form as Eq. (57) of |3q , this author approximated density peaks by density extrema. Therefore, our coefficients agree 
with his expressions only in the limit v 3> 1, in which b v — > v/cr® and b^ — > 0. As we will see shortly, the product 
of £> S pk(fc, zq) factors can be used to define spatial bias parameters at all orders. For peaks of significance v, the first 
order biasing is equivalent to the Fourier space multiplication by b\{k, zq) = 6 sp k(fcj zo), i.e. 

bi(k, z ) = bio + b 01 k 2 where &i = b v , b i = ■ ( 10 ) 

We emphasize that this result is exact: there are no higher powers such as k 4 , etc. In by, i and j count the number 
of factors of b v and b^, respectively [our notation should not be confounded with that of 87i]. In Scc. ^IIIBl we will 
demonstrate that the bio are the bias parameters in the local bias model. Eq. pJJl) defines the first order bias for peaks 
of height v. Notice that, in real space, \>\{k, Zq) = bio — boid 2 is a differential operator acting on fields and correlation 
functions. Hence, the first order average peak overabundance can also be rewritten <5 p k(^, x, zq) = (bi<5,s)(x, 20). 
Using the peak bias ^ , it is straightforward to show that the real space cross- and auto-power spectrum are 

P £s(^ R s, k, z ) = bi(fc, z ) Ps(k, z ) W(kR s ) (II) 
P$(y,Rs,k) = b 2 1 (k,z )P s (k,z )W 2 (kR s ) . (12) 
The corresponding relations for the correlation functions are 

eSU^S.r,*)) = (bid 0) ) = i>io^ \Rs,r,zo) + boi^ 1 \Rs,r,z Q ) (13) 
$(y,Rs,r) = (b?^ 0) ) =blod°XRs,r, z ) + 2b w b 01 ^\r s , r,z ) + b 2 01 ^\R s , r,z ) . (14) 



Note that the cross-correlations with the linear density field <5(x, zo) depend explicitly on zo- As shown in [7J,l75(, these 
expressions agree with those obtained from a rather lengthy derivation based on the peak constraint, which involves 
joint probability distributions of the density field and its derivatives. It is worth noticing that, while expressions (|12[) 
and (|14p are only valid at leading order, the cross-correlation functions (TlT|) and (TT3"]) are exact to all orders. 

We emphasize that the biasing ([5]) is a mean bias relation that does not contain any information about stochasticity. 
Due to the discrete nature of density peaks however, one can expect that the average peak overabundance <5 p k(x) in 
a cell centered at x generally be a random function of the underlying matter density (and its derivatives) in some 
neighborhood of that point. In fact, while the bias is deterministic in Fourier space, it is generally stochastic and 
scale-dependent in configuration space [75| . 
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C. Velocities 



In what follows, we will be interested in the gravitational evolution of the correlation of initial density peaks for 
which the velocity field also matters. As can be seen from, e.g., the average bias relation ([S]), peaks locally move with 
the dark matter (since the gradient of the density vanishes at the position of a peak). However, the three-dimensional 
velocity dispersion of peaks, er 2 pk = (Vp k ), is smaller than the mass velocity dispersion a?_ 1 [H, |H[, 

«4 k = ^ 1 (l- 7 ?) I (15) 

because large scale flows arc more likely to be directed towards peaks than to be oriented randomly. As recognized in 
[75j . the fc-dependence of the first order peak bias b spk (k) leads to a fc-dependence of the peak velocity statistics even 
though the peaks move with the dark matter flows. Taking the divergence of the peak velocity Eq. ([5]) and Fourier 
transforming, we find 

6 pk (R s , k, z ) = (l - ^ fc 2 ) W(kR s ) 9(k, zq) = b vpk (k) 6 S (K z ) , (16) 

where 6 = V • v is the velocity divergence. This defines the linear velocity bias factor 6 V pk(fc) for peaks of significance 
v and curvature u. Note that it does not depend on v, u nor on redshift and, for the highest peaks, it remains scale 
dependent even though the spatial bias bi(k, zq) has no fc-dependence. Nonetheless, for notational consistency, we 
define 

b V pk(fc) = & vp k(fc) = &vpk(fc) (17) 
as being the velocity bias of peaks of height v. 



D. Smoothing scale and peak height 

To illustrate the key predictions of the peak formalism, we will present results for the two-point correlation of 
density peaks in a ACDM cosmology with h = 0.7, f2 m = 0.279, fib = 0.0462, n s = 0.96 and normalization erg = 0.81 
consistent with the latest constraints from the CMB The sound horizon at recombination is ro f=s 105 h~ 1 Mpc. 

The peak height v and the filtering radius Rs could in principle be treated as two independent variables. However, 
in order to make as much connection with dark matter halos (and, to a lesser extent, galaxies) as possible, we assume 
that density maxima with height v = 5 sc (zq) / ao(Rs) identified in the smoothed density field linearly extrapolated at 
zq are related to dark matter halos of mass Ms oc R s collapsing at redshift zq, where S sc (zq) is the critical density 
for collapse in the spherical model [9(3, [Ml and we use a Gaussian filter to relate smoothing scale to mass. A more 
realistic treatment should include non-spherical collapse since the maxima of a Gaussian density field are inherently 
triaxial [9214951 ] . In the background cosmology we assume, the linear critical density for spherical collapse at zq = 0.3 is 
S sc rj 1.681. The Gaussian smoothing scale at which v = 1 is Rs,, ~ 0.8 /i~ 1 Mpc, which corresponds to a characteristic 
mass scale M St « 6.2 x 10 11 M Q /h. 

While there is a direct correspondence between massive halos in the evolved density field and the largest maxima of 
the initial density field, the extent to which galaxy-sized halos trace the initial density maxima is unclear. Therefore, 
we will only consider mass scales Ms significantly larger than the characteristic mass for clustering, Ms t , for which 
the peak model is expected to work best. For sake of illustration, we will present results for v = 2 (2a) and 
v = 3 (3a) density peaks. At redshift zq = 0.3, this corresponds to a filtering length Rs = 2.9 /i _1 Mpc and 
Rs = 5.3 /i _1 Mpc or, equivalently, a mass scale Ms = 3.0 x 10 13 M Q /h and 1.6 x 10 14 M /h. To help set scales 
in the discussion which follows, the associated values of (71, co/ui, bio, boi) are (0.65,3.7,0.8,21.1 h 2 Mpc -3 ) and 
(0.68,6.2,3.5,72.4 h 2 Mpc -3 ), respectively (note that bias factors here are Lagrangian ones). The three-dimensional 
velocity dispersion of these peaks is a 2 _ 1 (1 — Jq): for our two smoothing scales, this corresponds to (7.75 /i _1 Mpc) 2 
and (7.03 /i _1 Mpc) 2 (recall that our velocities are in units of aHf as 58 kms^ 1 /iMpc -1 at z = 0.3, so dispersions 
have dimensions of (length) 2 ). 



III. CORRELATION OF INITIAL DENSITY PEAKS AT SECOND ORDER 



A. The general formula 

Correlations of density maxima can be evaluated using the Kac-Rice formula [9(| [9?J ■ In this approach, r)i (x) is 
Taylor-expanded around the position x p k of a local maximum. The number density of peaks of height v' at position x 
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in the smoothed density field 6s reads as (we drop the subscript S in the right-hand side for notational convenience) 



n pk ( Z / / ,i? s ,x)=^|detC(x)|^ 3 )[77(x)]^[A 3 (x)]5Kx)-^ 



where 



1± 

(72 



(18) 



(19) 



is the characteristic radius of a peak. Note that Eq. (TTS)) is independent of the redshift zq for peaks at fixed v. The 
three-dimensional Dirac delta 6^ 3 \rj) ensures that all extrema are included. The product of the theta function 0(A3), 
where A3 is the lowest eigenvalue of the shear tensor Qj, and the Dirac delta 6{v — u') further restrict the set to 
density maxima with specific height v' . The 2-point correlation function for maxima of a given significance separated 
by a distance r = |r| = 1x2 — xi| thus is 



1 + £ P k(^, Rs,r) 



ipkjv, Rs, xi)n pk (i/, R s , x 2 )) 
"pk(^ R s) 



where n p k(i / , Rs) is the differential average number density of peaks of height v on filtering scale Rs |36 

1 



n v ]Ay, Rs) 



- 2 /2 G (i)( 7l)7lI/ ) 



(20) 



(21) 



{2ir) 2 R 3 

Note that it does not depend on zo (or, equivalently, on the amplitude of density fluctuations) at fixed v. The 

function Gq (71,71^) is defined in Eq. (|A60j) . For the 2a and 3cr density peaks considered here, the mean abundance 
is n p k = 9.0 x 10 -5 and 4.8 x 10~ 6 h 3 Mpc~ 3 , respectively. While the calculation of Eq. (|2T)|) at first order in the mass 
correlation and its derivatives is rather straightforward |7J| (this is Ea lM)) . at second order it is quite involved. The 
main steps are detailed in Appendix [A] Fortunately, most of the terms nicely combine together, and the final result 
can be recast into the compact form 



^Rs,r) = (bUi 0) ) 



-4(^ 1/2) bnd 1/2) 



2oj 
3 
2aJ 



10 



.(2)n2 
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( 2 h 2 l 



■»(a)/ 



(3 



(i)\ 2 



Uh 2 



3(f. 



(3/2)n2 



(3/2)n2 



(22) 



In the right-hand side of Eq. (f22j) . all the correlations are function of Rs, r and zq. More precisely, the first line 
contains terms involving first and second order peak bias parameters bi and bu, the second line has a ^-dependence 
through the function 1 + (2/5)9 Q lnGo (71, 7if)|a=i (which is displayed in Figj9|), and the last two terms depend on 
the separation r (and Rs) only. Note that this expression exhibits not only terms quadratic in bias parameters but, 

unlike standard local bias (Eulcrian or Lagrangian), also terms linear in them. These terms involve derivatives 

of the linear mass correlation £q that vanish at zero lag. They arise because the peak correlation depends also on 
the statistical properties of r\i and Qj . 

In analogy with bi(fc, zq), the action of the second order peak bias bu is defined as the Fourier space multiplication 
by 



bii(gi,92, z o) = & sp k(gi, zo)fr sp k(<72,zo)- (1 - 7x 



■bix(qi+qi) + b 02 q(q J 2 , (23) 



where q\ and qi are wavemodes and the coefficients 620, bu and bo2 describing the peak bias at second order are 



&2o(^, Rs,Zo) = K„ 
611(1/, Rs,z ) = b u(i + 
bo2(v,Rs,z Q ) = b cc - 



1 



1 



^(l-7 2 ) 
7? 

^(1-7?) 
1 

* 2 2 (l-7i 2 ) 



v — 2^\vu + j(u 



2 „2 



1 



C0O2 



(l-7i 2 ) 2 (!-7i 2 ) 
(l + 7 2 ) z/u - 71 [v 2 + M 2 ] 



7i 



(i - 7?) 

2 — 271 vu + jfiy 2 



(1-7?) 



(i-7?r 



(l~7i 2 ) 



(24) 
(25) 
(26) 
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FIG. 1: Left panel : Lagrangian bias coefficients characterizing the second order peak bias bn(gi, 92), Eqs. (|24[> ~ H26[) . as 
a function of peak height for a filtering radius Rs = 2.9 /i _1 Mpc or, equivalently, a mass scale M s = 3 x 10 13 M Q /h. The 
shape parameter is 71 w 0.65. For the 2<j peaks considered in subsequent illustrations, 620 is negative, 620 ~ —1.2. Right 
panel : The second and fourth root b^{ 2 and b^* define a characteristic scale below which the scale dependence of bn is large. 
In the limit v — s> 00, bo2 becomes negative and converges towards — ovT 2 (l — 7i)~ 1 , whereas bn asymptotes to the constant 
(7i/o"i) 2 (l — 7?)" 1 - Note that, in contrast to b\{ 2 and bf/ 2 4 that have units of length, 620 is dimensionless. 



See Fig. [T] for the relative size of these contributions as a function of peak height. As shown in Appendix fJA] 
b vv , byQ and arise upon averaging the products b u b v , b v b^ and b^b^ over the peak curvature. In this respect, 
= Gn \ii,1xv)/Gq (71 ,71^) is the nth moment of the peak curvature at a given significance v. For sake of 
completeness, bjj acts on the functions (t - ) and C|" 2 '( r ) according to 

-1 />oo />oo 

(dr^nd^) = 4^1 / d 9i y dg 2 g 2(ni+ V" 2+ % 2 i(9i,92^o^^ (27) 

When l\ = £2 = 0, the real space counterpart of bii(gi, qi, zq) is readily obtained by making the replacement q 2 — s- — d 2 
(which reflects the fact that ^q™'' (i?s, r, zq) is a solution to the Hclmholtz equation (d 2 + fc 2 )£(r) = 0). Eq. (j2"2")l is 
the main result of this Section. We note that [72j also computed second order corrections to the peak correlation 
£ p k(^, Rs, r ) for which, however, they did not provide any explicit expression. 

Before illustrating the impact of the second order terms on the correlation of initial density maxima, we remark 
that, although the calculation of the peak correlation at third order is very involved, the contribution proportional to 
(£o°^) 3 can be derived relatively easily. We find 



- 3b v 

uuu 



<7o 2 (l-7?) 

where the third order coefficient b,,,,„ is defined as 



m^U^)\ (28) 



- u 3 - 2> lx v 2 u + 2,-ijvu 2 - jfu 3 

b vvv {v,R s ,z ) = — -3 . (29) 

CT o (1 ~ Ti) 

Thus, up to third order, the peak correlation may be cast into the form 

fa(v,Rs,r) « bU { 0) + ibl (^ 0) ) 2 + ^ao^o ') 3 + additional terms , (30) 

where the missing terms, while of the same order as the ones we display, have a more complicated structure (see 
Eg. (|2"21 for second order contributions). In the limit v ^> 1, the scale-independent pieces bio, &20 and 630 to the 
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bias asymptote to the values bio — > v/oq, 620 — > (^/fo) 2 and &30 — > (^/co) 3 obtained in the high level excursion set 
approximation [35| . We will see shortly that these bias factors are indeed equal to the peak-background split biases 
derived from the average peak abundance Eq.([21"|). 

In Fig[TJ the second order Lagrangian biases 620, &11 and bo2 are shown as a function of the peak height for the 
mass scale Ms = 3.3 x 10 13 M Q //i (left panel). The second and fourth root b\{ 2 and bj^ 4 , respectively, define a 
characteristic comoving scale below which the corresponding scale-dependent terms bug 2 and bo2? 4 are large. In the 
limit v — > 00, the scale-independent piece b2o increasingly dominates whereas bn and bo2 tend towards the constant 
value (7i/cti) 2 (1 — 7 2 ) -1 and —a^ 2 ^ — 7i) _1 - For a realistic threshold height v < 4 however, the scale-dependent 
contributions cannot be neglected since |bn| and/or |b 2| are typically much larger than | &20 1 - Although the exact 
value of the second order biases somewhat changes with the mass scale Ms, their overall behavior varies little over 
the range Ms ~ 10 12 — 10 14 M Q //i as 71 weakly depends on Rs- Therefore, our conclusions hold regardless the exact 
amount of smoothing. It should also be noted that, if one wishes to associate these bias factors to halos of mass 
Ms = 3.3 x 10 13 M Q //i, then the variation with v is in fact a variation with rcdshift. 

The correlation £ p k(is, Rs, r) is shown in Figj2]for the 2a and 3cr initial density peaks collapsing at rcdshift zo = 0.3. 
The solid (green) curve represents the first order term b 2 £o°' (Ea.(fI3)l) while the long-dashed-dotted curve is the 
full second order correlation (Eq. (j2"2")) . We have also plotted the second order contributions quadratic in bn, linear 
in bn and independent of bn separately. They are shown as the short-dashed-dotted, short-dashed and long-dashed 
curve, respectively. Notice that (b2o, bn, bo2) = ( — 1.2,23,363) and (7.8,285,3927) for the low and high threshold, 
respectively. For the 3<r peaks, the term linear in bn is negative over the range of distances considered and, thus, 
appears as a dotted line. In fact, the piece linear in bn is the only negative contribution at small separations but 
it vanishes at zero lag. Since the true peak correlation rapidly converges to -1 for r < Ri (as shown by [7fl |. the 
small scale behavior of £ p k is dominated by an exponential term exp(— i? 2 /r 2 )), small scale exclusion should manifest 
itself in higher-order terms, but it is beyond the scope of this paper to calculate them. We can also observe that the 
correlation of 2a peaks is negative on scales r ~ 5 — 10 /j _1 Mpc. However, this is likely an artifact of truncating the 
expansion at second order in the correlations 

Figure [3] focuses on the baryon acoustic oscillation. At separation r > 80 /i _1 Mpc, second order corrections are 
negligibly small, so that £ p k(^, Rs, r) is given by Eq.(fT4"|) with an accuracy better than 1%. For comparison, we also 
plot the first order correlation bf £<s arising in a local biasing scheme with same value of bio- It is important to 
note that £5 is the correlation of the unsmoothed, linear mass density field (in practice we use Rs = 0.1 h~ Mpc). 
It is quite remarkable how the "sombrero" -shaped terms 2bioboi£o^ a nd boi£o restore the contrast of the baryonic 
feature otherwise smeared out by the large filtering (Recall that Rs = 2.9 and 5.3 /i _1 Mpc for the 2a and 3a peaks, 
respectively). The BAO in the peak correlations is even sharpened relative to the BAO in £5. A thorough discussion 
of this effect can be found in [7J]. In mV[ wc will see that, although most of the initial enhancement of the BAO 
contrast is smeared out by the gravitational motions of the peaks, some of it survives at the epoch of collapse. 



B. A peak-background split derivation of the peak bias factors 

The peak-background split [H, HH, H3, |98[ is a heuristic argument that furnishes another way to derive the large 
scale bias of density peaks. This approach is quite different from ours because it is based on number counts in 
configuration space and, thus, does not make reference to the bias in Fourier space. 

There are two ways in which the peak-background split is implemented. In the first [36T.l37l.l42l]. the iV-th order bias 
parameter 6jv( 1/ ) is related to the iV-th order derivative of the differential number density n(y) of virialized objects 
according to 



1 \ N _, ,-id N [n{v)] 
a {z )J dv ! 



b N (u > z )^(-- 7 - T ) niv)- 1 ^^, (31) 



with the important caveat that the mass function is universal (i.e. it depends on v solely). For density peaks, on 
setting n — n p k(^, i?s) and performing the derivatives with respect to v at fixed smoothing radius Rs, one obtains 


b N (u,R s ,Zo) = ( ^"^(^^' ^'^f 511 • (32) 

V a (Rs,zo)J ov N 

As noted by [75|], the first order peak-background split bias is 



bi(u,R s ,z ) = b 1Q (v,R s ,z ) 



(33) 
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FIG. 2: The correlation of initial density peaks at the second order is shown as the long dashed-dotted (magenta) curve for 2a 
(left panel) and 3<j (right panel) density peaks collapsing at redshift zq = 0.3 according to the spherical collapse prescription. 
For the Gaussian filter used in this paper, this corresponds to a mass scale Ms = 3 x 10 13 and 2 x 10 14 Mq//i, respectively. 
The individual contributions appearing in Eq. (|22[) are shown separately. Namely, the solid (cyan) curve is the first order 
contribution i>i^\ whereas the second order term quadratic in bn, linear in bn and independent of bn are shown as the short 
dashed-dotted, short-dashed and long-dashed curve, respectively. A dotted line indicates negative values. 



This means that the large scale, constant and deterministic bias factor returned by the peak-background split argument 
is exactly the same as in our approach, when we are on large enough scales that the fe-dependence associated with 
the boifc 2 term can be ignored. It turns out that Eg. ([55)1 generalizes as follows: higher order derivatives of the peak 
number density (|2ip with respect to v [which are reported in |42[ result in the large scale, fc-independent peak bias 
coefficients 



b n (v,R s ,z ) = b 2 o(v,Rs,zo)i 



nil 



(u, R S , Z ) = b 30 (f , Rs, Zq), 



etc. 



(34) 



However, derivatives of Eq. (f2lj) cannot produce the fc-dependent bias terms like brji, &n etc., which arise owing to 
the constraints imposed by derivatives of the mass density field. 

Therefore, we will now consider the second implementation of the peak background split [33l . l98j in which the 
dependence of the mass function on the overdensity of the background is derived explicitly. The ratio of this conditional 
mass function to the universal one is then expanded in powers of the background density. The bias factors are the 
coefficients of this expansion. We will demonstrate below that this is the correct approach to recover the scale or 
/c-dependcncc of the peak bias parameters. 

The key quantity is the average number density of peaks identified on scale Rs as a function of the overdensity 5b 
defined on another smoothing scale Rb (we use the subscript B because we are mainly be interested in the regime in 
which the scale Rb of the background satisfies Rb 3> Rs)- This conditional peak number density is 



i p k(^, Rs\Sb, Rb) 



f4 0) (7i,7i£) exp[-(y - ev B ) 2 /2{\ - e 2 )] 



(27r) 3 / 2 i?? 



(35) 



where 



vb 



5b_ 
cob 



(u\u,u B ) = 71 v = 7i f 



to = e = 



7i 



1 



Ox 



1-r 



(uv B ) = 71 er, 



Var(u|z/, vb) 



(fc 2 



- 2 

7i 



oTx/c? 



s 



O"0x/^, 2 



05 

2 (i 
1 



(36) 
(37) 
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with a n s = cr n (Rs, zq), a n B = <J n {R B , Zo), and we have defined 

1 



2 _ 



2^ 2 Jo 



dk k 2(n+1) P s (k) W(kR s )W{kR B ) 



(38) 



[see equation E5 of[3Q|. Here, the x denotes the splitting of smoothing scales, i.e. one filter is of size Rs, the other of 
size Rb- We have deliberately written r as an average of k 2 to emphasize that we naively expect it to give rise to the 
k 2 dependence of peak bias. This will eventually be proven correct. In addition, note that tvs = (8b/o~qs) ( (7 ox/ <t ob)- 
In what follows, it will be convenient to define £ 2 B = cr ox/ (T oB- When Rb 3> Rs, this ratio is of order unity (there 
is a form factor that depends on the shape of the smoothing filter). 

Notice that the integral of Eq. (|5"5j) over all 6b gives the unconditional number density n p k(^, -Rs) of Eq. ([2"TT) . The 
peak background split expands the ratio n p k (v, Rs\8 B , Rb) / ?ipk (v, Rs) hi powers oi 8b- This ratio is then interpreted 
as representing the average overabundance of peaks in regions which have mass overdensity 8b although, strictly 
speaking, it is a statement about cells of overdensity 8b that have a peak at their center. Therefore, it is not a 
statement about randomly placed cells, even though, as we discuss below, it is often treated as such. 

If we set r — > and e — > then the coefficient of the term of order 8b gives bio, that of order gives S B gives i>20, 
etc. Note that in this limit, 71 — > 71 and 71 v — > 71^(1 — Sb/Ss), so that expanding Eq. ([33|) in powers of 8b will be 
the same as differentiating Eq. (j2"Tj) with respect to 8s- These derivatives result in the large scale, fc-independent, peak 
bias factors we have been denoting as bio, b2o etc. As noted above, however, these derivatives cannot produce the 
fc-dependent bias terms like boi, bn, etc. 

Setting e — > but keeping the r dependence means that 71 — > 71 but jifi — > 7iz/[l — (1 — r)(8 B /8s)^ 



x b\ ■ As a 



result, derivatives of G °^ with respect to 8b will introduce terms which depend on r; these are terms which could not 
have been obtained by differentiating the unconditional mass function. For example, to first order in 8b, 



G 0) (7i,7i^«G 0) (7i,7 H 



l-(l-r) 

1 - 7i °s 



so that Eq. (f3"5j) becomes 

n p k{v, Rs\8b,Rb) ~ n pk (v, R s ) 
Hence, in this limit, 



t , S b v2 n \ ( 7i )7i b ^2 

1 + v E^ B - (1 - r) g 7- S x 

o-qs 1 - 7i °s 



(39) 



(40) 



(<5 p k|<5i 



_ n p k(v, Rs\5 b ,Rb) 1 V 2 
= ; ^ — ; 1 ~ 



n pk (^, i?s) 



(70S 



1-7? 



(5s 1 - 7i 



°"ox <5b 


^ - 7iu> 


+ 1 x 5B 


- 7i^\ 


^o(-Rb) CT 5 


U-7i 2 V 


^OB °2S 


Vl-7? J 



UX r c 1 " J-X t, jr 

OlO OS H 2~ 001 



(41) 



'OB 



'OB 



where bio and boi were defined in Eq. ([TUl) . Therefore, the cross correlation between the average overdensity of peaks 
defined on scale Rs and the mass overdensity on scale Rb is 



(<5 p k^s) (Rb) = o"ox bio + CT ix boi on large scales i? B > Rs- 



(42) 



Note that this final expression is the Fourier transform of (bio + boik 2 )Ps(k, zo)W(kRs)W(kRB)- Thus, we have 
shown explicitly that our implementation of the peak background split argument has produced the same linear, scale 
dependent bias as the Fourier space argument. 

At this point, one may be worried there is an inconsistency in the above expansion since e — > necessarily implies 
r — > 0. For example, for Gaussian smoot hing of a powerlaw spectrum P(k) oc k n , we find r ~ 2(Rs / Rb) 2 and 
e ~ (2Rs / Rb) 3+ "^ 2 at sufficiently large Rb [36|. For a realistic spectral index —3 < n < 3, e and, obviously, r always 
vanish in the limit Rb — ^ 00. In fact, there is no problem here since one can, at least formally, treat e and r as 
two independent, small parameters in addition to 8b- Hence, setting r — > at fixed value of 8b and e for instance 
corresponds to retaining the change in the small scale density Ss while neglecting the change in the curvature u 
induced by the background wave 8b- 

The higher order bias factors can be derived in an analogous way. Each term of order <5 B will include terms of 
order r m with m < N . Terms proportional to r give boi, bn, etc.; terms proportional to r 2 give bo2, 612, etc. Thus, 
our analysis provides a simple way of determining all the additional fc-dependent higher-order bias terms which arise 
in the peaks model. The fc-dependent polynomials associated to these bias factors are determined from symmetry 



11 




80 100 120 140 80 100 120 140 

r / h~'Mpc r / h _1 Mpc 



FIG. 3: A comparison between the initial unsmoothed density correlation £(r, zo) (black, dotted-dashed) and the initial peak 
correlation £ p k(*A Rs, r) (red, solid) around the BAO. To obtain the peak correlation, the density field was smoothed with a 
Gaussian filter on mass scale Ms — 3 x 10 13 (left panel) and 2 x 10 14 Mq//i (right panel). The dotted-long dashed, short-dashed 
and long-dashed curves represent the individual contributions bio£o°'> 2biot>oi£o 1 ' ) an< ^ ^oi^o^ to the nrs t or der peak correlation 
(Eq |14|l . A nonzero boi restores, and even amplifies the acoustic peak otherwise smeared out upon filtering the mass density 
field. The dotted curve indicates the second order correction to the peak correlation. Results are shown for the CDM transfer 
function considered in this paper. 



considerations. For instance, 612 multiplies the polynomial q\q\ + + q\q\- As a general rule, the bias factor by 
is associated to the monomonial symmetric function m(i 1; ... (Qi , ■ • • , if+j)- 

Our peak-background split derivation is very interesting for the following reasons. Firstly, all the analytic models of 
halo/galaxy clustering assume that the (scale-independent) bias coefficients multiplying powers of the mass correlation 

function £q (Eq. |3T)|) are the peak- background split biases (Eq. l3"Tj) . In the peak model, this equivalence can be derived 
from first principles. Secondly, the peak-background split holds even though the mean number density fipw{v, Rs) is 
not universal (due to its Rs dependence). For a spherical collapse prescription v = <5 S c(-2o)/cos, this implies that it 
is the abundance of virialized objects of mass Ms as a function of collapse redshift zq which is related to the bias 
parameters, and not the mass function - the abundance as a function of mass Ms - at fixed redshift zo- Thirdly, 
when correctly implemented, the second of our prescriptions yields the correct scale-dependent bias factors, despite 
the non-universality which comes from the fact that the background density is correlated with the curvature. And 
finally, the peak-background split approximation has been shown to provide a good description of large scale peak 
bias in simulations [42l. |99| although, recently, deviations at the 5-10% level have been reported [iil. Il0dj | . Since our 
expressions reproduce this limit, we have confidence that our approach will furnish a good approximation at the smaller 
scales where the bias parameters become scale-dependent (which we showed is reproduced by the peak-background 
split). 

Summarizing, we have shown how to determine the cross correlation between the peak point process and the 
smoothed, linear mass overdensity 5b order by order. In fact, because we have the full expression for the correlation 
between peaks and the surrounding density field, this cross-correlation can be computed exactly. The calculation 
simplifies by noting that the integral to be done is f(u) g{u\v, Vb) 9(v\vb) Kb g{ v B), where g is a Gaussian variate. 
This can be rewritten as f(u) g(i>) g(u\v) vb 9{vb\v, u). Assuming that vb ranges over all [—00,00] simplifies the 
integrals to compute, and one finds 

a pk s(Rb) = °ox &io + CT ix boi at all scales. (43) 

This expression, which is exact to all orders, coincides with the average mass profile around an initial density peak 
[3ol l75l ]. This average density profile has been shown to provide an accurate description of the cross-correlation 
between the Stokes Q and U polarization parameters and the temperature peaks of (two-dimensional) Gaussian CMB 
maps (HI. 
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C. Relation to the local bias model 



The local bias model is commonly used to model the clustering of dark matter halos and galaxies. In its Lagrangian 
formulation, this model assumes that the tracer overdensity 6^ of halos of mass Ms is a local deterministic function 
of the linear mass density field smoothed on scale Rb- This function may be written as the Taylor series [e.g.. l39j 

5 h (v,M s ,R B ,x) = b +J2 b N(») J ^r, (44) 

AT=1 

where the dependence on Ms arises through the peak-background split bias parameters b]\r(v) solely. The value of 
bo is set by requiring (5h(v, Ms, Rb-, x)) = 0. The choice of Rb is arbitrary as long as the constraint Rr ^> Rs is 
satisfied. The large scale clustering is then independent of the smoothing on scales kR B "C 1 (see [H, [55[ ; see also 
[5l| for a discussion of smoothing in the context of local Eulerian biasing) . Therefore, the cross correlation between 
the tracers overabundance and the mass overdensity in spheres of radius Rb is 

al s (R B ) = (5 B (x)5 h (v,M s ,R B ,x)) = £ b N (u) <M*)*fl(*)) w L + h ni0 * B \ a \ B . (45) 



A comparison with Eq. (l43l) reveals that, even on the largest scales, the cross correlation between the peaks and the 
mass distributions, a 2 k Sb ~ bja^ x , differs from the local bias expression, which involves all the bias factors with odd 
N values. This can be traced to the fact that the peaks expression is for cells centered on peaks, and only the mass 
field is smoothed on scale Rb-, whereas the local bias expansion is for randomly placed cells, where both the tracer 
field and the mass have been smoothed on scale Rb ■ 

For similar reasons, the auto-correlation function, at second order in the local bias model, is 

W,M 8 ,R B ,r) « 6? $\R B ,r,z ) + \b 2 u [$\R B ,r,z )] 2 , (46) 

where we have ignored extra terms involving powers of a^ B which arise owing to the continuous nature of the bias 
relation (|44]l. 

This demonstrates that, for large separations r 3> Rb, the peak correlation function £ p k(^> Rs, r ) is consistent with 
that of the local bias model although, for peaks, the bias factors are fc-dependent. It is pretty clear from Eqs. pp| and 
(|30|) that the local bias scheme is a special case of the more general peak model which, on large scales, approaches a 
local, deterministic, scale independent relation only in the high peak limit v 3> 1. Notwithstanding this, the previous 
analysis shows that 



i v v{v » 1, R s , r) « b]^\R s , r, z ) + ± [b 2 n + JJ^^^ d 2 ^ \^ ] (R S , r, z )\ + other terms . (47) 

on smaller scales where the fc-dependence of peak bias matters. As can be seen, the second order peak bias (the term 
in curly brackets) remains scale-dependent even for the highest peaks. Therefore, if the peaks model is correct, then 
we shall expect deviations from the local bias model on mildly nonlinear scales k ~ 0.01 — 0.1 /iMpc -1 , even for the 
most prominent density maxima. 
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IV. GRAVITATIONAL EVOLUTION OF THE PEAK CORRELATION FUNCTION 



Thus far, we have explored the scale dependence of bias in the 2-point correlation of local maxima of the primordial 
density field. However, nonlinear collapse and pairwise motions induced by gravitational instabilities will distort the 
initial correlation. Since the precise calculation of the dynamical evolution of £ p k(^ 7 Rs, r ) is rather involved, we 
will assume that the initial density peaks arc test particles which flow locally with the dark matter according to the 
Zcl'dovich approximation (i.e. peaks move along straight lines). We will calculate the correlation of their positions 
as a function of redshift and show that we recover a velocity dampi ng fa ctor and a mode-coupling power similar to 
that found in (Eulerian) renormalized perturbation theory [RPT, see Il0l| . 



A. The peak correlation function in the Zel'dovich approximation 

In the Lagrangian approach, the Eulerian comoving position and proper velocity of a density peak can generally be 
expressed as a mapping 

Xpk(z) = q P k + S(q pk ,z) , v pk (z) = a(z)S(q pk ,z) , (48) 
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where q p k is the initial peak position, S(q, z) is the displacement field, a is the scale factor and the peak velocity is 
in standard units (i.e. not scaled by aHf). A dot denotes a derivative with respect to cosmic time. We will assume 
that the local maxima are test particles that do not in terac t with each other. Therefore, at the first order, the peak 
position is described by the Zeldovich approximation jl02l |. in which the displacement factorizes into a time and a 
spatial component, 

S = -£>(»V$(q), S = -/3(z)V$(q) . (49) 

Here, D(z) is the growth factor of linear mass density perturbations and $(q) is the perturbation potential linearly 
extrapolated to present time. Explicitly, $(q) = </>(q, z)/AirGp m {z)a 2 D(z) where 0(q, z) is the Newtonian gravita- 
tional potential and p m {z) is the average matter density. Finally, /3(z) = HDf, where H(z) is the Hubble constant, 
is proportional to the l ogarithmic derivative / = dh\D / dh\a. Note that f(z) scales as £l m (z) - 6 for a wide range of 
CDM cosmologies (To^ |. 

Let us now consider an ensemble of realizations of initial peak positions. The correlation function £ pk (i/, R$,r, z) 
(we will henceforth omit the dependence on v and R$ for brevity) is related to the zeroth moment of the joint 
probability P 2 (vi, v 2 ; r, z|pk) to have a pair of peaks separated by a distance r and with normalized velocities vi and 
V2 . Following [l04| , we write 

nl k [l + t P k(r,z)] = y'd 3 v 1 d 3 V2P 2 (vi,v 2 ;r,z|pk) . (50) 

As we will see shortly, even though the probability P 2 ( vi , v 2 ; r , z | pk) depends upon the distance r = r • f and the unit 
direction vector f , the peak correlation depends only on r. When the peak motions are governed by the Zel'dovich 
approximation Eq. (|49|) . we can easily relate P 2 (vi, v 2 ; r, -z|pk) to the joint probability distribution at the initial redshift 
Zi > 1, 



P2(vi,v 2 ;r,z\pk) = P 2 (v 1 ,v 2 ;r - <r v Av 12 , Zi\pk) = Jd 3 r' <5 (3) (r' - r + a v Av 12 ) P 2 (v 1 , v 2 ; r', ^|pk) 



(51) 



where a v = <r_ 1 (z) is the 3-dimcnsional rms velocity variance of the matter, and Au 12 is the velocity difference v 2 — v\. 
Note that in this equation and those that follow, velocities have been scaled by aHfa-x, i.e., v = v/cr_i. Eq. ([5T|) is 
a cons equence of Liouvillc's theorem, which states that the phase space density of peaks is conserved [so, as shown 
in ll04L one can easily obtain a differential equation describing the evolution of the n-particle distribution functions]. 
It is especially useful because we know how to calculate 2-point distributions subjects to the peak constraint in the 
Gaussian initial conditions, 

Pa^t^r'^lpk) = Jd 6 C 1 d 6 C 2 n pk {^ 1 )n pk { X ' 2 )P 2 (w u w 2 ;r / , Zi ) . (52) 

Here, P 2 is a joint-probability for the 13-dimensional vector w = (v^, rji, v, £a) of variables at position and x.' 2 , 
and ripk(x) is the Klimontovitch density Eq. (fT5]) . Expressing the Dirac delta as the Fourier transform of a uniform 
distribution, we find 

n 2 pk [l+U(r,z)} (53) 
= Jd 3 r'e^ r -^ y'd 6 Cid 6 C 2 n pk (x' 1 )n pk (x^)|y'd 3 Wl d 3 ^ 2 P 2 (w 1 ,w 2 ;r',z l )e 4 ' T " k ' Aui2 | 

= / (£^ J d3r ' eik ' (r_r,) J d6 Cid 6 C2 n pk (x' 1 )n pk (x^)P 2 (y 1 , y 2 ; r', z t ) ^Jd 3 v 1 d i v 2 P 2 (v u v 2 \y u y 2 ; r', z t ) e ^ Av A 

where the vector y corresponds to (r/j, v 1 (a)- The second equality follows from Bayes' theorem. To integrate over the 
velocities, we use the identity 



Jd N yy n ■■■y in P(y) e a ' y = H) n -£- ' ' ' ^ ex P (-^ jtsj + ijts ) ' 



(54) 



which follows from the relation 

1 r r . / \ 1 1 

(55) 
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Here, P(y) is a iV-dimcnsional Gaussian multivariate of covariance matrix £ = (yy^) and centered at y = H. Eg . ([54]) 
is a very useful relation since it allows us to circumvent the inversion of the covariance matrix. On inserting the above 
identity into Eq. (|53|) . the peak correlation function may now be formulated as 



rr 



P k 



[1 + £ P k(r, z)} = J Jd 3 r' ^ Jd 6 Cid e C2 n pk (x' 1 )n pk (x^)P 2 (y 1 , y 2 ; r', z t ) cxp f-^EJ + iJ f S 



(56) 

where J = cr v (k, — k). The task of computing the rcdshift evolution of the peak correlation function boils down to the 
evaluation of the 6-dimcnsional covariance matrix S and mean vector 3. 



B. A simple illustration: the mass correlation function 



To understand the physical meaning of nonlinear corrections as well as emphasize the relation with, e.g., RPT, 
it is instructive to consider first the unbiased case, i.e. the evolution of the matter correlation function. Therefore, 
there is no peak constraint, so the two-particle probability density Piivi, r, z)d 3 V\d 3 Vi is simply the probability 
to find a pair of dark matter particles separated by a distance r and with velocities v\ and v%, respectively. As a 
consequence, £ is the covariance of matter velocity components, i.e. £y = (viVj), and S = because there is no net 
mean streaming between two randomly selected locations. After some simplifications, the mass correlation function 
evolved with the Zel'dovich ansatz takes the simple form 



1 + U(r,z) 



d 3 k 

(27T) 5 



e lk r e 3 fe 



d 3 r'e jkr exp 



Co 



(-i) 



-i) 



a v A-i) 

9 ?2 



(k.f) 



(57) 



This agrees with Eqs.(16) and (17) of |l04j provided that his <f>(r) corresponds to our £q (r), so that A</>(0) = — (T 2 _ 1 . 

It should be noted that tr_i and Q ^ are evaluated at redshift Zi ^ 1. Hence, the ratio a v /a i is equal to D(z) / 'D(zi). 
After some manipulations, the last exponent can be rc-cxpresscd as [in the notation of Il05| 
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On taking the Fourier transform of Eq. (|57p . we recover the well-known expression for the nonlinear mass power 
spectrum in the Zel'dovich approximation, 



P m {k,z) = e 
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(60) 
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where qi... n = qi H hq« and the kernels F n [1061 . Il07| are symmetric, homogeneous functions o f the wave vectors 

qi,...,q„ that describe the nonlinear evolution of the density field in the Zel'dovich approximation [l08| . 



K(qi,- 



1 (k • qi) (k • q„ 



qi 



q n 



(61) 



They are nearly independent of Q m and Oa- Note that, in the exact dynamics, F 2 (qi,q 2 ) = 5/7+ l/2(q±/q2 + 
<72/<7i)qi • Q2 + 2/7(qi • c^) 2 - In the Zel'dovich approximation however, the coefficients are all equal to 1/2, reflecting 
the fact that momentum is only conserved at first order. 

In Eq. (|60[) . the sum repres ents all the mode-coupling corrections, whereas the decaying exponential pre- factor 
corresponds to the propagator |l05| . which for Gaussian initial conditions des cribe s the (imperfect) correlation between 
the nonlinear and linear density field. Hence, the Lagrangian formulation of |l04| furnishes an easy way to obtain the 
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rcsummed RPT propagator [see also Il05l 1 103 . for a similar Lagrangian description]. At the second order, the mass 
correlation is 



1 e^e-? k2 °-^(p s (k,z) + -^ /d 3 qi /d 3 q 2 F 2 ( qi ,q 2 ) ' P s {q u z)P s (q 2 , z) <5^(k - qi - q 2 ; 



(2tt) 3 



(2tt) 3 



(62) 

The second order kernel scales as F 2 (qi,q 2 ) cx k 2 i n the (sq ueezed) limit where the total momentum k goes to zero 
as a consequence of mass-momentum conservation 



C. Including the peak constraint 



The computation of Eq. (|56[) is rather involved when the peak constraint is taken into account. Here, we will 
explain the basic result and show that it generalizes previous formulae based on the local bias scheme. Details of the 
calculation can be found in Appendix $B] where it is shown, among others, that it is quite convenient to work with 
the Fourier transforms of £ and H. We eventually arrive at 



£ P k(r, z) 
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r/ 3 k r i 2 

^3 G 2 (k, z) [bf (ft, z)\ P Ss (k, zo) e lkr + £ MC (fc, z) 



Ps s (k, Zi) + PMc(k, z) 



(63) 



The last equality follows upon making the replacement z,; — > zq (everywhere but in the exponential pre- factor). Here, 
zq is some fiducial rcdshift which we take to be the collapse redshift (because results are usually normalized to low 
redshift quantities). The function 



G 2 (k,z) =Dl(z)e-i k V 



(64) 



where D + (z) = D(z) / 'D(zq) is the linear growth rate normalized t o its value at the epoch of collapse, is a damping 
term induced by velocity diffusion around the mean displacement Il04|. Namel y, G(k, z) is the peaks propagator in 
the Zel'dovich approximation, and is analogous to that introduced in [10li [l05j for the matter evolution (the latter 
depends on u_i). In the exponent, the factor of 1/3 reflects the fact that the dynamics governing £ p k(?", z) is effectively 
one-dimensional (because, on average, only the streaming along the separation vector of a peak pair matters). The 
first order Eulerian and Lagrangian peak biases bf(k,z) and bi(fc, zo), both defined with respect to Ps s (k, zq), are 
then related according to 



bf(k, z) = b V p k (fc) + D+\z) bi(fc, z ) . 



(65) 



In the limit z — > oo, we recover bj(k, zo)Ps s (k, zq) at the first order. The presence of b vp k(fc) reflects the fact that peaks 
stream towards (or move apart from) each other in high (low) density environments, but this effect is fc-dependent 
owing to the statistical velocity bias. Still, the fc-dcpcndcnce of b vp k(fc) is such that the li near Eule rian peak bias bf 
is scale- independent in the limit ft <C 1, in agreement with the "local bias theorem" [55l 1 112L Ill3| |. On writing the 
Eulerian linear peak bias as bf = bf Q + b^k 2 , Eq. ([65|) becomes 



bf (z) = l + D-\z)b 10 (z ), 



b^(z) = D7\z)boi(zo)--%. 



(66) 



The first relation is the usual formula for the Eulerian, linear scale- independent bias [33|. It shows that, unsurprisingly, 
the bias of the peak distribution eventually relaxes to unity 114l4ll6| |. The second relation implies that bgi approaches 
the negative, i?s _ dependent constant — a 2 /a 2 as the gravitational instability grows. A scale dependence in the linear 
peak bias bf thus persists in the long term if the linear velocities are statistically biased (otherwise b^ — > as z — > oo). 

Thinking of the first order peak statistics as arising from the continuous bias relation Eq. (J5j) furnishes a straightfor- 
ward derivation of Eq. ([66|) . Namely, the linear continuity equation for the mass density field reads 5(t) = —V • v(t) = 
6{t) whereas, for density peaks, S p k(t) = —V • Vpk(f) = P k(t) = bvpkW^^- Recall that ^pk^ is an operator which, in 
Fourier space, is a multiplication by (1 — a 2 s /a\ s k 2 )W{kRs). The solution is 
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(67) 
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As a result, the leading-order contribution to the cross correlation between the average overdensity of peaks defined 
on scale i?g and the mass density field smoothed on an arbitrary scale Rb satisfies 

(68) 



(S pk (t Q )S B (t )) = \ D + (t)6 B (t)[(p+ \t) - l)b vpk W5(t) + S pk (t) 

= D~\t) \{D~\t) - 1) (s B (t)b vpk S s (t)) + (6 B (t)6 pk (t)) 

= D-\t) 



{D-\t)-l) (a 2 x (t)-^a 2 x (t)) + b 10 (t)a 2 x (t) + b i(t)a 2 x (t) 
\ a is J 



On writing the left-hand side as bio(£o)CT 2 x (to) + &oi(io) <T ix (*o) an d isolating the terms in cr 2 x and a\ x , we arrive at 



bxo(t o y ox (t ) = D+\t) (D+\t) - l)a 2 x (t) + bio(t)a 0x a 2 x (t) 



b 01 (t ) ( j 2 x (t ) = D- 1 (t) 



(D^(t)-l)^a 2 x (t) + b 01 (t)a 2 lx (t) 



IS 



(69) 
(70) 



Using the fact that the spectral moments of the linear mass density field scale as cr nX (to) = D + 1 (t)a nx (t), we can 
rewrite the above as 



bio(t) = 1 + D^(t)(bio(t ) - l), b i(t) = -^ + D+ 1 (t)(b i(to) + ^f-) , 
v ' a \s V a is / 



(71) 



which is precisely Eq. (|66|) . On linear scales, although the peak bias is deterministic in Fourier space, it is generally 
stochastic and scale dependent in configuration space [75| . In principle, it would be possible to solve also for the cross- 
correlation coefficient between the peaks and the smoo thed mass density field Sb by considering the time evolution 
of the peak rms variance (6 2 k (t)) for instance [see, e.g.. ll!5| . 

With the peak constraint, the one- loop contribution to the mode-coupling power Pyic(k,z) in the Zel'dovich ap- 
proximation eventually becomes 
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4(2tt)3 a\ s o\ s 
xPs s (q 1 )P Ss (q 2 )6^(k-q 1 -q 2 ) . 

where all spectral moments a n and power spectra Ps s (q) are evaluated at zq. The second order Eulerian bias 
bjj(qi, q 2 , z) is a symmetric function of q\ and q 2 , 

bfi(qi,q 2 ,z) = -7 r 2(qi,q 2 ) + ^D^(z) .Fi(<li)&i(?2, ^o) + ^lM^qi, z ) + ^D+ 2 (z)b~ u (qi, q 2 , z ) . (73) 

Here, bu(qi, q 2 , zq) represents the second order Lagrangian peak bias Eq. (|23p and, in analogy with standard PT, 
we have introduced the kernel JF n which characterize the nth-order evolution of the peak correlation function in the 
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FIG. 4: The evolved power spectrum P p ^(k, z) (dashed curve) for the 2a peaks as predicted by Eg. (|63|l at the redshift of collapse 
z — zo — 0.3 (left panel) and at z = 1 (right panel). The first order term and the mode-coupling contribution are shown as 
solid curves together with the unsmoothed linear mass power spectrum Ps(k,zo). At low wavenumber, the scale-dependent 
contribution to the mode-coupling power, Pmc(^) — Puic(k = 0), scales as fc 2 since there is no mass- momentum conservation. 
The Poisson expectation l/n p k (not shown on this figure) is « 10 4 . 



Zel'dovich approximation, 

Ai(qi) • • • On) = -7 bvpk 91) • • ■ bvpk(q'n) ■ (74) 

n\ qi q n 

J- n is identical to the standard PT kernels except for the velocity bias b V pk(<z)- As can be seen, the second order mode- 
coupling power is simply obtained from Eq.([2"2")l upon a replacement bn — > bn m Ea.((25]l. a Fourier transformation 
and a multiplication by the diffusion damping pre-factor exp(— (l/3)fc 2 cr^ pk ). Therefore, the mode-coupling induced 
by gravity is Eq. ([72")l minus the second order terms in Eq. ([2"2")l . We also note the plus sign in the second term of ([72)1 , 

which follows from the fact that the Fourier transform of £1 (r) is iq(r ■ q). Finally, we simply Fourier transform 
PMc(k, z) to obtain the mode-coupling contribution in configuration space. 

In Eq. ([56)l . each power of J brings a factor of Ti(q). At the second order, there is only one contribution proportional 
to J 4 ~ J~2 but, contrary to the Eulerian PT expression, it does not involve the first order bias. In fact, taking the 
local bias limit in which bi = &/, bn = bu and ignoring the exponential damping pre-factor and a possible statistical 
velocity bias (i.e., G = 1 and b vp k = 1), we do not recover the familiar PT expression [e.g.. |4(1] 

2 

P 5s (q 1 )P Ss (q2)S {3 Hk-c ll -q 2 ) , (75) 

where bf^ arc Eulerian peak-background split bias factors, but rather (omitting sub-leading powers in the growth 
factors for simplicity) 

2 

P Ss ( gi )P Ss (q 2 ) 8^ (k - qi - q 2 ) . (76) 

The difference is not surprising. In the first expression, local bias is applied to the evolved density field and, thus, 
bias coefficients enter as multiplicative factors in each term of the perturbation series, including the mode-coupling 
contribution. In our approach, however, peaks are identified in the initial conditions and then evolved gravitationally 
such that, in the long term, the Lagrangian bias factors drop to zero and the unbiased case is reproduced. In fact, 
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FIG. 5: Matter and peak correlation functions at the redshift of collapse z = zo as predicted by Egs. ()62p and (|63p . respectively. 
In configuration space, the mode-coupling contribution £mc is always positive left to the BAO and negative right to the BAO, 
so it introduces a shift in the acoustic scale. The magnitude of £mc depends strongly on the bias parameters. 



our Eq. (f76|) agrees with the local Lagrangian bias expression of |l!7l | , 

2 

P Ss ( qi )P Ss (q 2 )S^(k- qi - q2 ) , 

(77) 

except for a pre-factor (1 + b ) multiplying the second order kernel F 2 (qi,q2) (in the peak model, this pre- factor is 
unity since the zeroth order bias is bo = 0). Therefore, in the local bias limit considered here, the one- loop contribution 
to the mode-coupling can be thought of as arising from the (nonlocal) relation J!?} 

1 + «5 pk (x, z) = [1 + <5 pk (q)] [1 + <J(x, z)\ . (78) 

Here, 5 p k(x, z) and <5(x, z) are the Eulerian peak and matter overdensity, 5 p k(q) = bi5(q) + (l/2)6//<5(q) 2 is the 
Lagrangian peak overabundance and x = q + S(q, z) is the mapping from Lagrangian to Eulerian coordinates. 
We expect also our Eg. ([76]) to agree with the local Lagrangian bias expressions of |55[ if we include higher order 
perturbative corrections to the peak displacement. We defer the calculation of the mode-coupling at second order 
Lagrangian perturbation theory (2LPT) to a future work. 
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D. Mode-coupling power: shot- noise and shift in the baryonic acoustic scale 



To help visualize the results of the previous Section, Fig.0]shows the evolved power spectrum P p k(&, z) of 2a density 
peaks as computed from Eq. (|63p at redshift z = zq = 0.3 and z = 1. The peak power spectrum is split into two 
distinct parts. The first piece is the linear mass power spectrum P5 smoothed with the filter W(Rs, k) 2 , amplified by 
a scale-dependent pre-factor (b vp k(fc) + D^ 1 {z)b\{k, zq)) 2 and damped with a diffusion kernel. The second piece is a 
sum of the mode-coupling power present in the initial conditions (induced by the peak biasing) and of that generated 
during the nonlinear evolution (induced by gravity). 

At small wavenumber, the mode-coupling power arising from the peak biasing dominates and, in the limit k <C 1, 
contribute a pure white noise term whose amplitude is independent of redshift. Namely, the low-fc white noise tail 
is generated by the peak biasing and not by gravity. While it can be shown that the redistribution of the matter 
caus ed b y the nonlinear interactions cannot build a white noise tail Ps{k) cx fc° in the mass power spectrum [see, 
e.g.. 1 103l | . it is interesting that this holds also for tracers of the mass distribution (such as density peaks) which do 
not conserve momentum. In Eq. (|72j) . the redshift independence of Pmc(^ — 5 z ) is ensured by the fact that all the 
-D(z)-dependent terms come at least with a factor of J>,-i.2 which vanishes in the limit k — > 0. The amplitude of 
PMc(k = 0), however, strongly depends upon the peak height. A quick calculation shows that, as k — > 0, 
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FIG. 6: Redshift evolution of the correlation of 2a and 3a peaks collapsing at za = 0.3 as predicted by Eq. (|63l) . The curves 
from bottom to top represent £ p k(^, Rs, r, z) at redshift z — oo, 5, 2, 1, 0.5 (dotted curves) and z = zq (solid curve). Only the 
correlation at the collapse epoch (z = zo) can be measured in real data. For comparison, the dashed curves show the correlation 
at z = oo and z — zq in a local bias approximation (see text). 



All three integrals converge owing to the filtering of the mass density field. However, the first and third terms 
are always positive whereas the second term can be negative. For the 2a peaks, these are 3.2, -1.4 and 1.9 x 10 3 , 
respectively. As the threshold height is raised, the first term increasingly dominates (we find 10.6, -2.1 and 1.0 x 10 4 
for the 3(7 peaks). Therefore, at second order in the expansion, the shot-noise, which we define as the sum of Eq. (|79[) 
and the Poisson expectation l/n p k, is super-Poisson for v > 2. This predicti on seems to be at odds with recent lines 
of evidence suggesting that the shot-noise of massive halos is sub-Poisson [52l . lll8l[ll9| . However, we caution that this 
super-Poisson behavior may be an artifact of truncating the computation of the peak power spectrum at second order. 

(n) 

We expect that, as higher powers of Q are included in the description of the peak correlation, small scale exclusion 
will increase and eventually make the large scale shot-noise correction sub-Poisson [seeHH, for a rough estimate of the 
effect] . 

Before continuing, we note that any local nonlinear bia sing intro duces constant power at small wavenumber in 
addition to the conventional 1/n shot-noise fill l49l, l5ll. [52l HT3l . 1 12dj . At second order, this white noise contribution 
is always positive as it is given by the first integral of Eq. (|T9")) with bu(q,q,z ) replaced by bu and S smoothed on 
some arbitrary scale R. Therefore, in contrast to the prediction of the peak model, the magnitude of this term is not 
well defined in local biasing schemes owing to the freedom at filtering the mass density field. 

Ignoring the k° tail, the next-to-leading term should scale as k 2 since there is no local conservation of momentum 
(which would otherwise enforce a k 4 behavior). One can easily check that this is indeed the case by noticing that, in 
the low- A: limit. 



h h h 2 h 2 



Fi(qi) = - Mj Fi(qa) = fi+— (l - 2 M 2 ) +0(k 3 ), F 2 (q 1; q 2 ) = -fx+- (l - 2 M 2 ) +0(k 4 ) . (80) 
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As k — s- 0, this implies Fi(qi) + i ? i(q2) ~ (^/<7i) 2 (l — , so the first two terms of scale as k 2 at low wavenumber. 

In configuration space, the mode-coupling contribution £,Mc{ r , z ) changes sign across the acoustic scale ro w 
105 h~ 1 Mj>c because the dominant term involv es th e derivative ^^(r) of the mass correlation which is positive 
(negative) to the left (right) of the acoustic peak |l2l| . This leads to a shift Ar of the inferred acoustic scale towards 
small scales fmUmj . Typically, a 1% shift in th e acoustic scale translates to a fairly substantial ~4% bias in the 



estimated dark energy equation of state [e.g.. \124 

As shown in FigfSJ the magnitude of ^Mc( r j z ) strongly depends on the first and second order bias parameters. For 
the 2cr peaks, the strength of the mode-coupling contribution relative to the linear term is slightly smaller than 
for the matter (~ 1 percent), whereas for the 3a peaks, it is much larger (~ 5 percent) due to the large, positive 
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FIG. 7: Ratio of the peak correlation £ p k(^, Rs, r, z) to the mass correlation £ m (r, z) at redshift z = oo, 5, 2, 1, 0.5 and zo (at 
separation 60 /i _1 Mpc, curves from bottom to top in the left panel). For convenience, the ratio is normalized to bio( z ) so that 
the results approach unity as r — > oo. The dashed curve shows the quantity at redshift of collapse zq = 0.3 in the local bias 
approximation. At z = zo, the bias of density peaks exhibits a ~5% scale dependence across the baryonic acoustic feature. For 
the 2<r peaks however, most of this scale dependence is caused by the initial amplification of the BAO contrast relative to that 
of the mass (see Fig[3} whereas, for 3<r peaks, this scale dependence is mainly generated by the gravity mode-coupling. The 
vertical line denotes the position of the linear BAO feature. 

first and second order bias factors. To estimate the "physical" shift of the acoustic peak generated by mode-coupling, 
we must correc t beforeh and for the "apparent" shift induced by the convolution with the propagator G 2 [we use the 
terminol ogy of 1 1 2 ll . Il22j . The latter is usually taken into account in the data analysis and cosmological parameter 
forecast |124lll25l |. Here, we will simply ignore the diffusion damping both in the linear and mode-coupling piece and, 
therefore, calculate the shift of the acoustic peak in the perfectly de-convolved correlation function 

£ p \>, z ) = [ bf ] 2 ® £<j 0) (r, z ) + £^ c (r, z ) . (81) 

It is important to note that, in the peak model, the "physical" shift is generated by second order mode-coupling and 
by scale dependence in the linear biases bi and b vp k- The relative contribution of the second source of shift, which 
is not present in the local bias approximation, turns out to be equally important at high and low peak height (This 
is somewhat counterintuitive since the scale dependence is more pronounced at low v). For the 2a and 3<r peaks, we 
find that the location of the acoustic peak in the linear mass correlation is "physically" shifted by Aro = —0.25 and 
-1.65 h~ 1 Mpc, respectively. The relative contribution of the first order scale dependence is in both cases ~ 20% only, 
mainly because bi docs not shift the Fourier phases. We believe these values should change somewhat at second order 
in the peak displacement (2LPT) since the cross term between the monopole and dipolc (of [t>n] 2 ), which generates 
most of the shift, will not remain the same. 

E. scale dependence across the acoustic peak: theoretical predictions and comparison with simulations 

We will now present results for the evolved peak correlation function and compare them to the auto-correlation of 
simulated dark matter halos. 

To begin, Fig. [5] shows the redshift evolution of the correlation £ p k(^, Rs, r , z ) of 2er and 3cr peaks from the initial 
conditions at z — oo (bottom dotted curve) until collapse at z — zq = 0.3 (top solid curve). The intermediate redshift 
values are z = 5, 2, 1 and 0.5. It should be noted that only the correlation at the collapse epoch can be measured in real 
data (assuming that the tracers are observed at the epoch their host dark matter halos collapse). For comparison, 
the bottom and top dashed curves represent the initial and final correlation for the local bias approximation in 
which (bi, b V pk) = (bi, 1) and the mode-coupling power is given by Eq. ([76]) . As the redshift decreases, gravitational 
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instability generates coherent motions which amplify the large scale amplitude of the peak correlation, and random 
motions which increasingly smear out the initial BAO feature. Although velocity diffusion due to large scale flows is 
less important for the peaks than for the locally biased tracers (owing to the fact that b vp k < 1), the final correlations 
are noticeably more similar than they were initially. Still, mild differences subsist at z = zq between the peak and 
local bias predictions, especially around the baryonic acoustic feature. 

In order to quantify these deviations, we take the square root of the ratio between the peak correlation £ p k (Eql63p 
and the mass correlation £ m ('Ea l62[l . Both arc consistently evolved at second order with the Zcl'dovich approximation. 
\/£pk/£m measures the scale dependence of peak bias as a function of separation. For illustration purposes, we 
normalize this ratio to bf (z) so that, on scales much larger than the acoustic scale (not shown in this figure), the 
normalized ratio rapidly converges to unity. Results are shown in FigJT] for the 2a and 3a peaks. The dashed curve 
represent the local bias prediction at the collapse redshift z = zq. Across the baryonic acoustic feature, the peak bias 
exhibits a residual scale dependence of ~5% amplitude. For the 2a peaks, this scale dependence arises principally 
from the initial amplification of the BAO contrast (see Fig(3J) whereas, for the 3a peaks, it is mostly induced by the 
gravity mode-coupling. In stark contrast to the peak model, the local bias approximation predicts negligible scale 
dependence for 2a tracers (there is a sharp upturn at r ~ 130 h~ 1 Mpc due to the fact that the zero crossings of £ p k and 
£ m are different). At 3a however, the discrepancy between both models is relatively smaller because the contribution 
of the mode-coupling generated during gravitational evolution, which is weakly sensitive to the bias factors boi, bn 
and bo2 in the limit v ^> 1, dominates the scale dependence of bias. The mode-coupling also contributes to suppress 
the peak bias by 2-3% at separations r ~ 60 — 80 h~ 1 Mpc [see also IrO] . 

Clearly, a very large simulated volume is required to search for similar scale dependences in the bias of the most 
massive objects created by gravitational collapse. Hence, we will present measurements of the baryonic acoustic 
feature in the auto-correlation of halos extracted from a single realization of 2048 3 particles in a cubical box of side 
7.68 /i _1 Gpc. The simulated volume thus is more than sixteen times the Hubble volume. Halos were subsequently 
found using a friends-of-friends algorithm with a linking length of 0.2 times the mean intcrparticle distance, leading 
to a final catalog of more than 15 million halos of mass larger than w 7 x 10 13 M©//i (with 20 parti cles or mo re). 
This N-body run and the associated halo catalog is a part of the MICE simulations project [see El EH, EH for 
an exhaustive description of the runs and Ihttp: / / www.ice.cat / mice for publicly available data] . A more detailed 



clustering analysis of this simulation will be presented elsewhere |145j . 

The filled symbols in FigJ5] show the measured bias (defined as the squared root ratio of their auto-correlation 
function to that of the dark matter field) for halos with significance larger than v t =2 and 3. Here, we exceptionally 
adopted a tophat filter to define the peak significance d sc (z)/ao. Therefore, the v t threshold corresponds to a mass 
cut M t — 2 x 10 14 and 7 x 10 14 M@/h, respectively (these halos have at leas t 35 and 200 particles respectively). 
The correlation was computed by direct pair counting using the estimator of [l28j |. Measurements were done first 
dividing the full boxsizc of 452 h 3 Gpc -3 into 27 non-overlapping regions of equal volume. The bias shown is the 
average over these regions and the error bars those associated with the corresponding variance of the auto correlation 
measurements (we actually depict the error on the mean, i.e., <7b/\/27). 

From Fig (5J it is clear that, while the > 3a halo sample is too sparse to furnish useful information, the correlation 
of > 2a halos shows unambiguous evidence for a scale dependence of bias of ^5% relative magnitude around the 
baryonic acoustic feature. To the best o f our knowledge, this is the first time that such a scale dependence is reported 
from N-body simulations [see, however. Il29j . 

To compare predictions from the peak model with the measurements from the N-body simulation, we remark that 
the set of dark matter halos with significance larger than v t includes all halos of mass larger than some threshold 
value Mt- As an attempt to reproduce the number counts of halos above a mass cut Mt, one may want to smooth 
the mass density field with a range of filter radii Rs and identify peaks of height S sc (zq) / ao(Rs) with virializcd halos 
of mass M = Ms- We adopt a sharp threshold even though the actual selection function w{M\v, Rs), which gives 
the probability that a peak of height v in the linear density field smoothed on scale Rs forms a virialized halo of 
mass M, may be fairly smooth [see, e.g., [3(1 for more realistic selection functions]. Consequently, we approximate the 
correlation function of dark matter halos as 



oo yoo 



£b(r,z) = N j—t / dvn pk (v,R s (v)) £ p k(r,z) where Nb(v t ) = j dvn^(y,R s (v)) . (82) 



v t J v t 



Here, Rs(y) is the filtering radius at which S sc (z) / ao(Rs , z) = v. The fraction of mass f{y)dv (or mass function) in 
peaks of height v thus is 

f{v)du = (M s /p) n pk (v, R s {v))dv , (83) 

and M s /p = (27r) 3 / 2 i?| is the mass enclosed in the Gaussian filter. In our simulation, the measured cumulative 
number density for the > 2a and > 3a halos is Nh « 1.45 x 10 -5 and sa 5.10 x 10 -7 h 3 Mpc -3 , respectively. For 
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FIG. 8: Scale-dependence of bias around the baryonic acoustic feature of the halo auto-correlation. Filled symbols with error 
bars represent the measurement for > 2a and > 3a halos extracted from the MICE N-body simulation (see text) at z = 0. 
The corresponding mass cut are M > 1.2 x 10 14 and M > 7.0 x 10 14 Mq//i, respectively. The solid curve corresponds to a 
prediction based on the peak model with linear peak-background split bias bi similar to that of the simulated halo catalogs. 
For the low mass cut, the threshold height is vt = 2.15 whereas, for the high mass cut, it is v t = 3 (see Eq |82p . A vortical line 
denotes the position of the linear BAO feature. 

these threshold heights, Eq.(|8"2"]l predicts slightly larger densities, i.e., Ah ~ 1-77 x 1CT 5 and ps 8.40 x 1CT 7 . To 
further improve the agreement with the simulated mass function, one shall ensure that a peak of height S sc on a 
smoothing scale i?s is not embedded in a region of height <5 SC on any larger smoothing scale. Carefully accounting for 
clouds in clouds in the peak formalism is a nontrivial problem even th ough, in a first approximation, we may simply 
enforce that the height of the peak be less than <5 SC on scale Rs + dRs [130L Il3l| . Fig.l of J8(J demonstrates that this 
substantially improves the agreement with the measured halo abundances at v > 3 but underestimates the counts at 
v < 2. Furthermore, adding this extra constraint also modifies the peak bias factors and, possibly, the bn-indepcndcnt 
terms in Eq.(f22j). For these reasons, we hereafter stick to the naive albeit reasonably good approximation Eq. (l82l) . 
and defer a more detailed modeling of the halo mass function to a future work. 

The predicted halo correlation Eq. ([52"]) is shown in Fig. [5] as the solid curve. For the low mass cut (left panel), the 
theoretical prediction assumes in fact v t = 2.15, which furnishes a better match to the linear, scale-independent halo 
bias measured in the simulation. For the high mass cut, v t — 3 as in the halo sample. Overall, since we consider the 
clustering of high peaks at z = 0, the scale dependence induced by gravity mode-coupling is important. Consequently, 
all the theoretical curves exhibit a broad peak structure centered at r ~ 90 — 100 /i _1 Mpc. By contrast, Lagrangian 
peak biasing solely would induces a "sombrero" -like scale dependence with a maximum close to the acoustic scale 
(see Fig(7]). Focusing on the \ow-v halo sample, it is difficult to assess the quality of the fit without knowing the data 
covariances. Nevertheless, the measurements strongly suggest that our prediction based on the peak model is a very 
good approximation. An improved description of peak motions, collapse and mass function will somewhat modify the 
theoretical curve but, because the predicted abundances are close to the measured ones, we believe the change should 
not be dramatic. It will be interesting to search for a similar effect in the clustering of low mass and/or high rcdshift 
halos, for which the contribution of the peak biasing should be larger. 

V. CONCLUSION 

In standard approaches to clustering, biased tracers are treated as a continuous field. Their overabundance is 
commonly assumed to be a local function of the smoothed mass density field. Gravitational motions can be straight- 
forwardly included through a spherical collapse prescription or perturbation theory to predict galaxy and dark matter 
halo correlation functions. 

In the peak model where the tracers form a well-behaved point process, the peak constraint renders the calculation 
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almost prohibitive. For this reason, all analytic studies of peak correlations thus far [36|, |38|, l72|, |7J, |75| have obtained 
results which are strictly valid at leading order and in the (Gaussian) initial conditions solely . 

To make predictions which can be directly compared with observational data or the outcome of numerical experi- 
ments, we have calculated the correlation of density peaks in the Gaussian linear mass field and computed its redshift 
evolution, consistently within the Zel'dovich approximation. We have then used our model to explore the effect of 
peak biasing on the baryonic acoustic feature (BAO). 

In the first part of this work (Sec. 3111) 1 . we obtained a compact expression for the correlation of initial density 
maxima up to second order. We demonstrated that the fc-independent pieces of the peak bias factors - bio, e ^ c - ~~ 
are equal to the peak-background split biases derived from the peak number density n p ^{v 1 Rs), even though the later 
is not universal. Furthermore, we also showed that the peak-background split approach can be applied to derive the 
fc-dependent bias factors - boi, &11, etc. This is especially interesting as it illustrates how the peak-background split 
should be implemented in any other model where the abundance of tracers depends on variables other than the local 
mass overdensity. So far however, we have not found a straightforward derivation of the terms linear and independent 
of bn in the initial peak correlation. 

In the second part (Sec. ^IV|) . we computed the gravitational evolution of the peak correlation under the assumption 
that the peaks locally flow with the matter. For simplicity, we considered the Zel'dovich approximat ion i n which the 
initial density peaks execute a motion along straight line trajectories. Using the approach laid out in Il04ll we s howed 
that the evolved peak correlation function can be expressed, like in renormalized perturbation theory 1 1 Oil . 1 1 9j as the 
sum of a linear contribution dominant at large scales, and a nonlinear contribution which describes the generation of 
power at a given scale due to the coupling between two modes (at one-loop level in our calculation) at different scales. 
The linear term can be easily derived from a continuity equation argument for the average peak overabundance if one 
includes the peak velocity bias b vp k- If one ignores the peak constraint, then our results reduce to those predicted 
by previous local Lagrangian bias calculations. Thus, the peak model is a useful generalization of this model. In 
particular, local Lagrangian biasing provides a good approximation to our results only in the limit of large peak height 
and separation, and only if velocity bias is ignored. 

We applied our model to predict the shape and amplitude of the baryon acoustic oscillation in the correlation of 



peaks at the epoch of collapse. Most of the initially strong scale-dependent bias across the BAO reported in [74 1 
is washed out by velocity diffusion - which manifests itself as an exponential damping kernel. Still, for 2er peaks 
collapsing at z = 0.3, our model predicts a residual ~5% scale-dependent bias which should be quite asymmetric 
about the acoustic scale. This prediction stands in stark contrast to the negligible scale dependence predicted by 
the local bias approximation for the same peak height. For 3<r peaks, the mode-coupling power dominates so that 
the predictions of both models are quite similar: the bias exhibits a broad bump at distances r ~ 90 — 100 ft. _1 Mpc 
smaller than the acoustic scale. 

We then measured the clustering of massive halos extracted from a very large N-body simulation. For halos with 
mass M > 2 x 10 14 M Q //i - which corresponds roughly to a peak significance v > 2 - the large scale bias shows 
strong evidence for a ^5% scale dependence around the baryonic acoustic feature. We have tried to reproduce this 
measurement using our model combined with a simple prescription to account for the observed halo counts. Our 
prediction, which is a weighted contribution of the scale dependence induced by Lagrangian peak bias and gravity 
mode-coupling, is in good agreement with the data. However, the model slightly overestimates the cumulative halo 
abundance. 

Clearly, there are a number of important missing ingredients - such as the second order contribution to the peak 
displacement and a better treatment of nonlinear collapse which must be accounted to achieve more accurate 
predictions. In this regards, all the machinery developed for the local bias scheme can be applied to the peak model, 
with the caveat that the bias factors, including the linear one, are /c-dependent. Numerical simulation s wi ll be 
essential to calibrate some of the parameters - such as the diffusion scale appearing the propagator [see, e.g.. ll32| - or 
determine what is the actual shape of the filtering kernel. Future work should also include redshift-space distortions, 
non-Gaussian initial conditions and massive neutrino species, and investigate the effect of small scale exclusion and 
stochasticity on the peak power spectrum and correlation function. 
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Appendix A: Correlation of initial density peaks at the second order 

In this Appendix, we present the derivation of the correlation of initial density maxima Eq. (|20[) at the second order 
in the mass correlation £q and its derivatives. 



1. SVT (scalar- vector-tensor) decomposition of the covariance matrix 

Computing £pk(?") requires knowledge of the joint probability distribution P 2 (y(xi), y(x 2 ); r) where y T = 
Ca) is a ten-dimensional vector whose components £4, A = l,--- ,6 symbolize the independent entries 
ij = 11,22,33,12,13,23 of Qj- 

We shall proceed in a way analogous to cosmological perturbation theory [e.g. Il33l4l40l and references therein], 
and decompose the variables into irreducible components according to their transformation properties under spatial 
rotations. Although this decomposition is not adequate for imposing the peak constraint, it greatly simplifies the 
structure of the covariance matrix C(r) of the joint probability density P 2 [see, e.g., [72[. In this work, we choose the 
covariant helicity basis (e + ,f,e_) as reference frame, where 

ie^ - eg . , ie^ + eg 

e+ = — — = — , r = r/r. e = — — (Al) 

and eg, are orthonormal vectors in spherical coordinates (0, <p). The orthogonality relations between these vectors 
are e± • e± = r • f = 1 and e + ■ e = e± ■ r = 0, where the inner product between two vectors u and v is defined as 
u • v = UiUi = UiV 1 . An ovcrline will denote complex conjugation throughout this Section. The property u ■ v = v ■ u 
follows from our definition of the inner product. For the first derivative of the density field, the SVT decomposition is 

77 = r/ s > + rf^ = r/ (0) f + ? / (+1) e+ + . (A2) 

Here, rf°> = tj ris the contravariant component of an irrotational vector (spin-0), r /\r]( s > = 0, whereas -q^ 1 ^ = rj e± 
are the two independent contravariant components of the transverse vector rf v ^ (spin-1), ry> v ' r = 0. The correlation 
properties of j/ - 1 and r]^ 1 ' can be obtained by projecting out the scalar and vector parts of the correlation of the 
Cartesian components r\i which, because of statistical isotropy and symmetry in i,j, is of the form 

( 7/l (x 1 )ry J (x 2 )) = H||(r)^ + H^r)^ , (A3) 

where 

Pij = (I3 - f (8 f) tf = e +l e +j + e-ie-j (A4) 

is the projection operator onto the plane perpendicular to f , and I„ is the n x n identity matrix. The vectors e± will 
be called contravariant basis vectors. They satisfy e± • e± = e±ie±i = 0, which follows from the relations e± = e±i 
and e±i = e' ± between contravariant and covariant basis. The covariant components of rf^' thus are rj = r\ • e±, 
whereas rj^ = rj^ . The functions H\\ = - 2^ ] )/2,(jI and H± = +(,2 ] )/3crj arc radial and transverse 
correlations [see, e.g.. Il4l| . Using n^ 1 ^ = r\ie l ± and rf^ 1 ^ 1 = T]ie±, we find that the correlations of the spin-0 and 
spin-1 components read 

(v i i ) V i 2 ° ) )=H ll (r), fap 1 ^) = H ± (r), (r/PVa^) = - (A5) 

Here and henceforth, the subscripts "1" and "2" will denote variables evaluated at position xi and x 2 for shorthand 
convenience. Let us now consider the symmetric tensor Qj. Writing Qj as the sum of a traceless symmetric matrix 
Qj and its orthogonal comple men t, Qj = Cy — (w/3)l3 where u = — tr^ = —Q, a suitable parameterization of the SVT 
decomposition for Qj is [e.g.. Il42| 



cu - -jus* + ^c (s) + v s 1 + ^ )f ) + (A6) 
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The functions u = — tr£ = — Q and = are the longitudinal and transverse spin-0 modes, Q are the 



components of a spin-1 vector, • f = 0, and Qj' is a symmetric, traceless, transverse (spin-2) tensor, S^qT 



Qj f J =0. Explicit expressions for these functions are 

C (5) = |5"» (lOT = 1 ( 3 r'r m - <, im ) Q m (A7) 
Cf 5 = v^Om = V5 $ - hf l ) f m Cim (A8) 

C? = ^pClm = \j\ [P\PT ~ \PijP lm ) Clm ■ (A9) 

We have introduced factors of ^/l/3 and \/2/3 in the decomposition (|A6[) such that the zero-point moments of the 
spin-0, spin-1 and spin-2 tensor variables all equal 1/5 (see equations (|A15[) . (|A16[) and (|A17[) below). Note that 
P( = e + ie J + + e~ie : L and P IJ = e!|_e+ +e l _eL to ensure nilpote ncy o f the operator P. The 2-point correlation function 
of any isotropic, symmetric tensor field Qj is of the form [e.g., Il43| 

(Cii( x i)0m(x 2 )) = Z\{r) fifjfif m + Z 2 (r) (fif[6 jm + fif m Sji + fjfi6 

+ Z 3 {r) (fifjS lm + fif m 6ij) + Z A {r) 6ijS lm + Z 5 (r) (8u5 jm + S im Sji) . 

In the case of the tensor Qj = didjS considered here, we have Zi = Z 3 and Z± = Z^. To project out the transverse 
spin-0, spin-1 and spin-2 parts of the 4-rank correlation tensor eq. (|A10[) . we act with the scalar, vector and tensor 
projection operators S a b, V^ c , Tjfij on (^(xi)^ m (x 2 )) and obtain 

(C (s) (x!)C (S) (x2)) = Z 1 (r)+4Z 3 (r) + 3Z 5 (r) (All) 

(C (y) ( Xl ) Cf } (x 2 )) = SPtj [Z 3 (r) + Z B (r)] (A12) 

(<$P (*0 CZ ] (^)) = 3T ijln Z B (r) , (A13) 

where 

Zl + 4Zs + 3Zs _ £ (I f - - f & + ftf >) , Z 3 + * - £ (±#> - itf> - ±tf >) , (M4) 

Z5 -^|il5 €o + 2p 2 + 35^ 4 

The correlation for the spin-0 variable Q 1 - ^ thus is 

(Cfcf) =^i+4Z 3 + 3Z 5 . (A15) 

Taking the 2 independent components of C^* 1 and their complex conjugate as being Q^ ±1 "> = ■ e± = y/3 e!±PQj 
and C (±1) = C (V} ■ e± = v^f^ yields 

(Ci^Cp 3 ) = Hz 3 + z 5 ), (Cx (±1) cf x) ) = o . (A16) 

Similarly, we choose C (±2) = Cif ^±4 = v/3/2 4 &j and C (±2) = Cy^ei = ^/372 e* ± e£Cy for the two indepen- 
dent polarizations and their complex conjugate, respectively. The correlation properties of these variables arc 

(d ±2) C 2 ±2) ) = 3Z 5 , (d ±2) C 2 T2) ) = . (A17) 

As expected, these are the only spin-2 degrees of freedom since Qj e l + eL = 0. Therefore, the second rank tensor Qj is 

fully characterized by the set of variables {u, Q°\ Q ±:L \ C^ 2 "*}- Furthermore, all the correlations are real despite the 
fact that some of these variables are c omp lex. We note that, in the particular case f = z, our variables are directly 
related to the variables yf m defined in [72(, who work in the spherical basis 

iy — x „ iy + x\ 

V- z ' V) ■ (A18) 
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The transformation from the helicity to the spherical basis vectors is performed by the rotation operator D\0, 8, <j>) 
(whose matrix elements are the Wigner D-functions 13^^,(0, 9, </>)). For the traceless tensor Qj, the spin-0, spin-1 and 
spin- 2 components in the spherical basis simplify to 



C (0) = \ (2C33 - Cn - C22) , C (±1) = T\J\ (Ci3 ± < 23 ) , C (±2) = yf (C11 - C22 ± 2zCi 2 



(A19) 



The relationship between the two sets of variables thus is £^ m )(x) = ?/2m( x )/v / 5 (see eq. (22) of [zll]). For the gradient 
rj(x) of the density field, the correspondence is rj^ m \x) = i/i m (x)/^/3. 

The 20-dimensional covariance matrix C(r) = (yy^), where y = (yi,y2) and y^ is its conjugate transpose, de- 
scribes the correlations of the fields at position Xi and X2. For simplicity, we will assume in what follows that these 
are smoothed on the same mass scale. C(r) may be partitioned into four 10 x 10 block matrices, the zero-point con- 
tribution M in the top left and bottom right corners, and the cross-correlation matrix B(r) and its transpose in the 
bottom left and top right corners, respectively. In terms of the variables yt = (r)^ ,Vi,Ui, £^ , , Cj > Vi ? Ci 1 7 
Ci +2) ,Ci~ 2) )> M and B have the block diagonal decomposition M = diag(M< ), M^, M (2 \ M^) and B = 
diag(B(°),B( 1 ), BW,B( 2 ),B( 2 )). Explicitly, 



/ 1/3 

1 71 

71 1 

V 1/5 



1/5 j ' 5 



(A20) 



and 



B(°) = 



3^lSo 

1 ^(1/2) 

1 c(3/2) 

CT1CT2 '1 
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(3^3 

B (l) 
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(3/2) >, 



1 r (V2) 
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^sf } 
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Off. 
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(A21) 

(A22) 
(A23) 



The second and third entries of the first row of B^ - 1 are the correlations of 772°'* with V\ and u\, respectively. As 
expected, these correlations are negative if i>\ > or m > since, in this case, the line-of-sight derivative is 
preferentially directed towards xi. The determinant of C(r) reads as 



detC(r) = detM (0) det 



M (o) _ B (o)" 



(m<°>) 1 



B (o) 



n 

s=l,2 



detM^] det 



jyr( s ) _ B (s)T 



(A24) 



It is worth noticing that, although C(r) docs not depend on the direction f of the separation vector r, it is not equal 
to the angular average covariance matrix C(r) = (l/47r) Jdfli- C(r). The latter follows upon setting Q n ' = whenever 
j in equations (|A20[) and (|A21|) . Furthermore, whereas the 2-point probability distribution i^yii y2! r) associated 
to the correlation matrix C(r) cannot be easily expressed in closed form, the joint probability density Pa(yij Y2] r) of 
covariance C(r) may be exactly written as the product 



■Pa(yi,y2; r) = P 2 (i / i,wi, v 2 , u 2 \ r) P 2 {r] 1 ,T] 2 - : r) P 2 (Ci, C2 



r) 



(A25) 



This factorization property reflects the fact that, upon angle averaging, the I — {v%,Ui), I = 1 (r/j) and £ = 2 (Qj) 
representations of SO (3) decouple from each other. Consequently, it should be possible to cast the 2-point probability 
densities in terms of rotational invariants such as the scalar product and the matrix trace. After some manipulations, 
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we find the joint probability density for the £ = 1 and 1 = 2 variables is 



^2 Oh: »72 ; r ) 



1 15 6 
20 (2tt) e 



(l-Tf)- 3/2 exp 



^l-T^)- 5/2 exp 



3ril + i-ql - 6T 1 r] 1 ■ r] 2 
2(1-T?) 

15\ tr(Ci 2 ) + tr (Cf)-2T 2 tr (ftfr) 

tJ (r^Tj) 



where we have defined T n (r) = (r) /c 2 for sake of conciseness. The following relations 



(o) (o) . (+i)_(+i) . (-1)_(-1) 

vi m + m m + m % -vi-V2 



AO) AO) 

Si S2 



E (d + ^ +s) + d- s) cr s) ) = |tr( Cl c 2 



(A26) 
(A27) 

(A28) 
(A29) 



can be useful to derive equations (|A26[) and (|A27|) . When the peak constraint is enforced, P2(ViiV2i r ) reduces to a 
simple multiplicative factor. Finally, the joint probability density P2(vi, u\,V2,u%; r) for the I — degrees of freedom 
evaluates to 

ftKm,^,m;r) = ■ , (A30) 

(2tt)V(1 -7i) A P 

Here, $(fi, «i, 1^2, "2; f) is the quadratic form associated to the inverse covariance matrix C~^(r) = (P,R T ;R, P). 
The 2x2 matrix P in the top left and bottom right corners reads as 



1 fPn Pi2 



Ap V Pl2 P 2 2 




-71 



1 1=^ 

7iTi(T -7i ! Ti)-7iT 2 (To-T 1 ) 

1=^ 



, 7iTi(T -7i 2 Ti)-7iT 2 (T -T 1 ) 

-71 + ^ ^ 

, Tg-2 7l 2 T T 1 +7i 2 T 2 
1 1=^ 



whereas the matrix R in the top right and bottom left corners is 



R = — 



1 



T - 7l 2 Ti) P n - 7l (T a - TO P 12 (T - 7^) P 12 - 71 (T a - T x ) P 
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(1 - 7l 2 ) A P V (To - 7i 2 Ti) P12 - 7i (T a - T x ) P 22 (T a - 7?Ti) P 22 - 7i (T - T x ) P ia 
Notice that the determinant Ap = P11P22 — P12 asymptotes to 1 — in the limit r->oo. 



(A31) 



(A32) 



2. Series expansion of the joint probability density 



To calculate the correlation function of initial density peaks at the second order, we first separate the covariance 
matrix into C(r) = C(r) + SC(r). The angular average C(r) = (M,B T ;B, M) contains the zero-point moments and 
the cross-correlation entries with £q ' solely, whereas 8C(r) = (O10, <5B T ; SB, O10), where O10 is the 10 x 10 zero matrix, 
encodes the cross- correlations We have for instance B^ 2 ) = ^ 2) /(5ct|) and 5B^ = (10^ 2) +3& ) )/(35a%). Using 
the identity det(I + X) = 1 + trX + (1/2) [(trX) 2 — tr(X 2 )] + • ■ ■ , we expand the joint density Pa(yi, y 2 ; r ) in the small 
perturbation 5C(r) and arrive at 



■P2(yi,y 2 ;r) « P 2 (yi,y2;r) 



1 - -tr( CT^C 



-tr^C-^CC-^c) 



(A33) 



1 + ^C-^CCHy + i (ytC-^CC-V) * - ^ytC-^CC-^CCHy 



at second order in SC. Here, the product C 1 <5CC 1 is of order 0(£) in the correlation functions Q , whereas 
trvC-^C), t^C-^CC-^C) and C^OC-^CC" 1 are of order 0(£ 2 ). Expressing the matrices C and SC in terms 
of the auto- and cross-covariances yields 



f 2 (yi, y 2 ; r) « P 2 (yi, y 2 ; r)| 1 + y^r^B M'Vi + \ (yl^SB M'Vi 



(A34) 
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at order 0(£ 2 ), where 

Q =2 (m -1 BM -1 <$BM -1 ) +M- 1 5B t M- 1 5BM- 1 . (A35) 

Note that <5B is not symmetric, so one must distinguish between 5B and its transpose. 

Let us consider the term y2M _x #BM _x yi linear in the correlation functions. Owing to the block-diagonal na- 
ture of C, it is a sum of contributions from the spin-0, spin-1 and spin-2 degrees of freedom. While the matrix 
(M w ) _1 5B w (MW) _1 fors = 1,2 generally has non- vanishing elements, it is easy to check that (M(°)) -1 (5B( )(M( )) _1 
has zero entries for the elements ij = 22, 23, 32 and 33. Furthermore, imposing the constraint r\ 1 = r} 2 = implies 
that y^ )t (M (s) ) _x 5B (a) (MW) _x y^ ) , s = 0,1,2 contains only terms linear in Ci\ & and products of the form 
Ci ■ At this point, we shall remember that the principal axes of the tensors Ci = C( x i) an d C2 = C( x 2) ar( 3 
not necessarily aligned with those of the coordinate frame. Without loss of generality, we can write (1 = — RAR T , 
where R is an orthogonal matrix that contains the angular variables (e.g Euler angles) and A is the diagonal matrix 
consisting of the three ordered eigenvalues Ai > A2 > A3 of — Ci- The value of u- t = — tr£i is invariant under rotations 
of the principal axes, while £1 transforms in the same manner as the £ — 2 cigenfunctions of the (orbital) angular 
momentum operator, i.e. the spherical harmonics Y^l 2 {v). Namely, on inspecting a table of spherical harmonics in 
Cartesian coordinates, we can write 



Y " {i) = lH ( 3z2 - r2 )' r 2 Y 2 ±1 (v) = T^d£z(x±iy), r 2 y± 2 (r) = - J £ (x ± iy) 2 . (A36) 



A comparison with eq. (|A19p shows that £( m ) ~ • x /5/47rF 2 m (f ). Therefore, the variables C}- m ' must transform in 
accordance with 

C (m') (x) =^D^ m ,(^tf,V)C (m) (x) (A37) 

tn 

under rotations of the principal axis frame. Here, Z? 2 nm /((/?, ip) ar e quadrupole Wigner D-functions with the Eulcr 
angles (tp,i),ip) as argument, whereas and £( in ) are the components of £ in the original and final eigenvector 
frames, respectively. Therefore, averaging over distinct orientations of the principal axes gives (£("*)) = 0. Noticing 
that, at the zeroth order, the joint density ^(yi, y2! ^) factorizes into the product -Pi(yi)-Pi(y2) of one-point proba- 
bility densities Pi(yi) which do not depend upon R, we find the correlations Q , j 7^ do not contribute to the peak 
correlation at the first order. This is in agreement with the findings of [72l |74||. 

The second order terms in the right-hand side of Eq. (|A34[) will also yield products in the variables £•{ and Q m ^ 
whose angle average can be reduced using the orthogonality conditions 

dRD i\ im ^ D%r4<P> = ^1 W^m,*,™^ (A38) 

d ^ D l m ^^) D t 2 rn^^) = (-ir^^^W-rW-miK ' (A39) 



SO(3) 

SO(3) 



Clearly, cross-products of the form Ci^C^™ ^ vanish upon averaging over the principal axis frames because they involve 
distinct rotation operators. However, when the angle average is taken at a single position xi = X2, we obtain 

77117712 

- f(™i)7( m *) X X , , 



5 

mi 77*2 

— 5 - / tr 

{ft °™ l l rn 2 L1 



where (• • •) denotes the average over orientations and we have also omitted the arguments of the Wigner /^-functions 
for brevity. Similarly, it is easy to show that 

^ C W) C «)J) = ^(mi)£K)^ = JL (-l)™^_ m , m , tr (c 2 ) . (A41) 
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These relations can be used to integrate out the orientation of the two eigenframes in the series expansion of the joint 
probability density i-2(yij Yi\ r )- For example, let us consider the contribution (y^M^B M _1 yi)(y2M _1 <5BM _1 yi). 
After some algebra and with the aid of Eq.(22) of [zH, we can write 



(ylM^BM-ViXyjM-^BM-Vi) = 



5d 0) d 0) + 5E (cl +s) c 
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s=±2 



(A42) 



on enforcing the constraint r7 x = r/ 2 = 0. To average over the orientation of the tensors Ci and C2, we note that 
products of the form f ("O^O^"*)^"*) simp lif y to 



>(mi)T(rai)>(m2)7(m2) 
SI S2 Si S2 
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tr C? tr C 



(A43) 



After some further manipulation, this leads to the cancellation of the term (|A42I) . 

However, the matrix traces in Eq. (|A34[) do not vanish on integrating over the angular variables, 
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nor does the second order contribution 



2 P 2(yi,y2;r) 



yjM-^M-Vi) -(ylQyi + y 2 Qy 2 ) « ^ l(yi) ^ l(y2) (y^aBM-Vi) -(ylQyi + y 2 Qy 2 

(A45) 

To evaluate the latter, we set the first derivatives to zero and recast the term yj,M _1 (5B M _1 yi; whose explicit 
expression is enclosed inside the curly bracket in the right-hand side of Eq. (| A42|) , into the following compact form 

ytM-^BM-Vi^^WC^+AMCf+SoMC^Cf+SiM £ d s) C 2 S) +52M £ ^ , (A46) 



s=±l 



s=±2 



where /1, / 2 , 50, <?i and g 2 are functions of r and, possibly, also u% and Mi (the exact expressions can be read off from 
Eq |A42[) . Upon squaring Eq. (|A46[ ). taking the average over the principal axis frames and substituting the relations 
(IA40P and (|A41|) . we arrive at 



y|M _1 5BM _1 yi 
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/i 2 tr(C2 2 ) + /2 2 tr(Ci 2 ) + ( ^ ) (.9o 2 + 2.9? + 2g 2 2 ) tr(C 2 ) tr(C: 
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(A47) 

fc(r)tr(<?)tr(c|) , 
(A48) 



The variables 6^ and 6^ arc defined as 



\ (Vi- 71 M, 



Co V 1 - 7f 



, hi = — 



1 f Ui- l\Vi 



°2 V 1 - 7l 



(A49) 
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They characterize the large scale bias of density peaks of significance z/j and curvature u%. As we will see shortly, 
product of these two variables generate bias parameters beyond first order. Likewise, the angular average of the 
scalar- valued function (yjQyi + y|Qy 2 ) can eventually be expressed as 

(yjQyi + y 2 Qy 2 ) = 4, (b vl & ] + b cl £, 2) ) 2 + ^ (b ul ^ 1/2) + b a ef /2) ) 2 + g(r) tr(cf) + 1 o 2 , (Aso) 

where 

(A51) 

is a function of the separation r solely. 



3. Second order approximation to the peak correlation function 

At this point, we follow [36| and transform the eigenvalues of C( x i) an d C( x 2) to the new set of variables 
{ui,Vi, Wi,i = 1,2}. Here, Vi and Wi are shape parameters that characterize the asymmetry of the density pro- 
file in the neighborhood of density maxima. After some algebra, the 2-point correlation of density peaks at second 
order in Q^to can ^ e written as 



1 5 5 3 4 (1 -T?)" 3/2 (1-T2)' 



5/2 



n% (2tt)6 i?6 v /(l- 7 2)A P 
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(A52) 
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2^ k (27T) 

5 



K b (1 - 7? 
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(3^ 2 + w 2 ) (3« 2 + w 2 ) h(r) 



+ JL (3 V ? + w l _ 1) (ft^CD + b(2 ^ _ jL (fr^CV*) + ^(3/2)) _ I ^2 + ^ q{r) + 1 ^ 2 

In what follows, we will restrict ourselves to the cross-correlation Cpk('r) = ?pk(^i, f 2 , Rs, t) of peaks of height v\ and 
j/ 2 . Therefore, we must integrate over the peak curvatures m. The peak constraint implies that the integration at 
fixed Ui > is restricted to the interior of the triangle bounded by (vi,Wi) = (0,0), (u;/4, — Uj/4) and (m/2,Ui/2). 
Moreover, is the quadratic form that appears in the 1-point probability density Pi(yj), 



2$„, = ut + 

1 



+ 5 (3w? + wf) 



v? - w?) 



(A53) 



(A54) 



(A55) 



Here, the integration domain is < ip < 2ir, 0<$<ir,Q<ip<2ir and dR = (l/8ir 2 )dcos , dd(pdip is the normalized 
Haar measure (J dR = 1) on the group SO(3). There is no analytic, closed-form solution to the integ ral Ip(Ai, A 2 ), 
although it can still be expressed as a hyper-geometric series in the argument (3Ai and A2 [see, e.g.,[l44j, and references 
therein] . We emphasize that equation (|A52[) is better than an approximation based on a second order Taylor expansion 
of the probability density -P 2 (yi,y 2 ;r) since it retains the isotropic part at all orders. 
For convenience, we write the peak correlation up to second order as follows: 



F(ui,Vi,Wi) is the weight function defined as [3 

F(ui,Vi,Wi) = (ui-2wi) (ui + w t Y - 9vf u 4 ( 
and2" (3 (Ai,A 2 ), with /3(T 2 ) = (15/2)T 2 /(1 - Tf), is the integral 

Zfl(Ai,A 2 ) 



/ dR exp 


/3tr(AiRA 2 R T ) 
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(A56) 
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where is the first order piece, cq. (fT4|) . and £^ are distinct second order contributions depending on i) correlation 
functions only ii) the peak height v and iii) linear and 2nd order bias parameters. We will now detail each of these 
contributions. 

(211 ^ 

• The first terms, £p k , follows from expanding the determinant of the covariant matrix C(r) at the second order. 
We have 



'1-7? 



T 2 - 4 7l 2 T T 1 + 2 7l 2 T 2 + 2 7l 4 T 2 + 2 7l 2 T T 2 - 4 7l 2 T!T 2 + T 2 

2(l-7i 2 ) 2 



(A57) 



Including the contribution from the trace, Eq. (|A44[ ). and expanding the multiplicative factor (1 — T 2 ) 3 / 2 (l 
T|)~ 5 / 2 « 1 + (3/2)T 2 + (5/2)T 2 at second order yields 



T 2 - 47 2 ToT! + 2 7l 2 T 2 + 2 7 4 T 2 + 2 7l 2 T T 2 - 4 7l 2 T!T 2 + T 2 
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(A58) 



(21) 

Notice that Q k does not depend upon the peak height, though it depends on the filtering scale Rs at which the 
peaks arc identified. 

(22) 

• The second contribution, Q k , contains all the terms for which the ^-dependence cannot be expressed as a 
polynomial in the linear and 2nd order bias parameters (to be defined shortly). For subsequent use, we introduce the 
auxiliary function 
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(A59) 



and its integral over the nth power of the peak curvature u times the u-dcpendent part of the one-point probability 
distribution, 



oo e -(x-w) 2 /2(l--fj) 

dx x n fix, a) 

v / MH) 



(A60) 



These functions are very similar, albeit more general than those defined in Eqs (A15) and (A19) of [36|]. With the 
above, moments of the peak curvature can now conveniently be written as u n (y) = Gn (7i j 7i i/ )/Gq ( 7 i,7ii / ). 

Next, we collect all second order terms that features the product of binomials (3f 2 + w 2 ) ni (?>v 2 + u> 2 )™ 2 with 
ni + n 2 < 2. For instance, expanding the integrand of Ip about j3 = gives 
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2^ 
45 



(3vj + w 2 ) (3v 2 + w 2 ) 



5T 



{3v 2 + w 2 ) (3v 2 + w 2 ) 



Similarly, 
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On inserting these expansions into the integral over the asymmetry parameters and using eq. (|A59[) , we find 
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1 5 
peak height v 



FIG. 9: The function 1 + (2/5)9 Q lnG*^ ) (71, 7i^)| a =i for several values of the correlation strength 71 = 0.2, 0.4, 0.6 and 0. 



This function vanishes in the limit v — > 00 since the logarithmic derivative of Gq tends towards —5/2 



where d a f(v,i, 1) = d a f(iii, a)\ a= i. We are subtracting the zeroth order contribution from the left-hand side of 
(|A63|) . which would otherwise give unity upon an integration over the peak curvature. There are two additional 
terms proportional to (3vf + wf)q(r) and a term like (3v\ + w\ )(3uf + vj\)h(r). For these terms, integrating out the 
asymmetry parameters yields 
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Adding Eqs (|A63|) and (|A64[) . we can eventually express the second order contribution Q]I ( r ) to the peak correlation 



function as 
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(A65) 



The last equality follows from the well-known relation Eq. (|2~Tj) for the differential density of peaks of height v. The 
logarithmic derivative of G^ with respect to a must be evaluated numerically. Nevertheless, it is worth noticing 
that g[, q) (7i,u;) and d a G { Q a) (71 , lo) are sharply peaked around their maximum. For large values of w, the former 
asymptotes to Gq ~ a -5 / 2 ix> 3 . Hence, this implies <9 Q lnGQ Q ^ (71, w)\ a =i w —5/2 in the limit u> ^> 1. 

• The last contribution, is the sum of two parts: a term which arises from the exponential e~* and, thus, 

(n) . . . (n) 

involves the angle average correlations £q , and a second part which involves the correlations QJq- Up° n expanding 
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e * at the second order and integrating over the shape parameters Wi, Vi and the peak curvature Ui, we obtain (after 
much tedious algebra) 
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where the second order bias parameters £>„„i = £w(fi,7i), 6^; = &j,f(i/i, 7 i) and b^i = b^ (^,71) are constructed 
from products of the variables b v and 6^ defined in Eq. (|A49[) , 
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Here, the overline designates the average over the peak curvature. Note that the last line in the right-hand side of 



eq. ([A66p is the contribution from the correlations 



4. A compact expression 

The second order contribution t;^ = £^ + £^ + may be written down in compact form with the aid of the 
second order peak bias operator bm defined through the Fourier space relation 



hli (91,32) = b spki (qi)b spki (q 2 ) - (l - 7i) 



(A70) 



where 6 sp ki(<7) = b vi + b^q 2 and q\ and q 2 arc wavemodes. By definition, its action on the functions Q (r) and 
it\r) is 
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The second order terms can be rearranged so as to recast the 2-point correlation of peaks of height v\ and ^2 into the 
more compact form 
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Notice that the first term of £^ combines with the first curly bracket of to give the second term in the right- 
hand side of Eq. (|A72)> . whereas the last term in Eq. (|X72]) is the sum of (3/2) T? and (3/crf )(^ 1) ) 2 in ^\ We recover 

Eq. ([22| in the particular case v± = 1*2 = v. For sake of illustration, the function 1 + (2/5)d a laG < ^\'ji,jiv) is shown 
in Fig. [5] for several values of 71. Note that it decreases monotonically and vanishes in the limit v — > 00. 

In the local bias model, the iV-order bias parameters are related to the iVth-order derivative of the mass function 
n(y) through a peak-background split argument (see ijlll Bl for details). Setting v\ = v-2 = v for simplicity and 
collecting the second order terms proportional to (f^) 2 that are present in and we find their sum is 
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Here, bjj is the second order peak-background split bias. Even though we do not calculate the peak correlation at the 
third order, it is straightforward to compute the coefficient multiplying (^q ' 1 ) 3 . This term arises from e~* at order 
0(£ 3 ) in the correlation functions, and e~* at times \/(l — 7 2 /Ap at C>(£ 2 ). Adding these two contributions 

yields 
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where bm is the third order peak-background split bias and the coefficient b vvv = bf, is defined in Eq. (|29)) . This 
demonstrates that the peak-background split derivation of the scale-independent peak bias factors holds at least up 
to third order. We speculate that this remains true at higher order. 



Appendix B: Nonlinear evolution of the peak correlation function 



In this Appendix, we provide technical details regarding the nonlinear evolution of the correlation of initial density 
peaks under gravitational instabilities. Here again, we will consider two different populations of density peaks of 
height v\ and vi identified on the same filtering scale Rs- 



1. The joint velocity distribution 



Computing the redshift evolution of the peak correlation function requires knowledge of the conditional probability 
density P%{v 1, U2|yi) y%', r, z%) which, in light of the SVT decomposition described in the previous Section, is 

P2(^ 2 |yi,y2;rU)= n p * (^ m \4 m VHytV^) ■ (bi) 

m=0,±l 

Here and henceforth, we omit writing the dependence on the initial redshift for brevity. Furthermore, the velocities 
v are in unit of aHFa-i (so there are dimensionless). Clearly, the line-of-sight velocity v^- ' = v ■ r' is a spin-0 
component that couples only to the other spin-0 variables y- ^ = y(°)(xj) = (r?j° , Vi, Ui, Q° ), while the components 
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^(i 1 ) — v . e± f the transverse vector correlate only with ^ = y' ±1 - ) (xi) = ,Ci ) an d the corresponding 

conjugates at x 2 . Using Schur's identities, we can write the conditional covariance matrix and mean value s( m ) 
for the velocity components as 



y(m 



_ y(m)T ^Q(m)~j 1 y(r 
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Here, is the covariance of the spin-m velocity components v^ m \ is the cross-covariance between and 

the vector y( m ) = (yj m \ y^" 1 ^), and C'" 1 * 1 = (y( TO )y( m )t) is the covariance matrix of y( m ). Like the 20-dimensional 
matrix C(r), C*- m ) can also be partitioned into four block matrices, with the auto-covariance M^™) along the diagonal 
and the cross-covariance 3^ m ' and its transpose in the bottom left and top right corners, respectively. The matrices 
MM and BM are defined in Eqs. (|A20[) and (|A21[) . To calculate the mode-coupling power, it is quite convenient 
to work with the Fourier transform of £( m ) and S^™-* rather than their real space counterparts (this allows us to 
circumvent the addition of three angular momenta). To this purpose, we will Fourier transform the entries of W (m) , 
B( m ) and CK m ) so that X( m )(r) = (1/8tt 3 ) / d 3 qX' m '(q) e 4q r . For the spin-0 variables, we obtain 
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/ 



\ 2a\ai 



k (r • a? 

__*2_(r. q ) 
-Ji-(r . q ) 

3 (f • q) 3 - (f • q) 



CTqCTI 



2(T0CT2 



(f -q) 

q 2 /a <7 2 
3 (f • q) 2 - 1 



— (r • q) 

q 2 /<7 a 2 



2a 1 <7 2 
I 



q 

2ai 



3(r -q) 2 -! 



2<7o(X2 



<i 
■Irr 1 , 



3 (f • q) - (r • q) 



3(f -q) 2 -! 



3(f -q) 2 -l 
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and 



7o/3 



V(°) T (q) = 
where 0i X 3 = (0,0,0) and 



0i 



X3 



(f-q) 2 V T 



7o/3 0ix3 



W (0) (q) = 



1/3 

■ (r • q)' 



3(f -q) 2 -l 

- (r • q) 

'l/3 



1 



»3 



(f • q) , — (f • q 



a 2 



2cr 2 



3 (f • q) - (f • q) 



For the spin-1 variables, the cross-covariances are 

& (e± • q) (e± • q) 



V3iq 3 



(r • q) (e± • q) (e± • q) 



(TifT2 

3q' 



(r ■ q) (e± • q) (e± • q) 



(f • q) 2 (e± ■ q) (e± ■ q) 
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(B4) 
(B5) 

(B6) 
(B7) 




7o/3 
(e± • q) (e± • q) 



\/3iq 



o 

q) ( e ± 



q) (e± ■ q) 



W (±1) (q) = 



1/3 

(e± ■ q) (e± • q) 



(e± • q) (e± • q) 

7o/3 

(e± • q) (e± • q) 

1/3 



^(f-q) (e ± -q) (e± • q) 




(B8) 
(B9) 



We have omitted a factor of Ps{q) in all matrix elements for shorthand purposes. Furthermore, we have not substituted 
r • q by q^> and e± • q by ' to avoid heavy notation for their complex conjugates. We also note that entries 
involving a odd power of f change sign under the space reflection f — > — f . 

With the SVT decomposition introduced above, the integral over the peak velocities can be written 



kA " 12 =exp( — jtEJ+ijtS 



(BIO) 



where £ = diag(E^°- ) , E^ -1 ') and 3 = 'ES +1 \ H* -1 )) are the covariance and mean of the multivariate Gaus- 

sian; and J, are the 6-dimensional vector a v (k^°\ —k^ ', k^ +1 \ —k^ +1 \ k^ x \ — fc'" 1 )) and its conjugate transpose, 
respectively. We will now Taylor expand the right-hand side of Eq. (|B10[) around the zeroth order contribution to 
and E< m ). 
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2. Zeroth order: diffusion damping due to random velocities 

At the zeroth order, the covariance matrices of the spin-0 and spin-1 conditional velocity distributions reduce to 

£ (0) = £ (±1) ^(l-7o 2 )l 2 , (BH) 
whereas there is no net transverse velocity S at this order. Inserting this result in Eq. (|B10j) yields 

yd Ar yP 2 (v 1 ,v 2 |y 1 ,y 2 ;r')e^ k - Av - « e^* 2 ^ 1 "^) , (B12) 

where &l pk (z) = a v( z ){^ " 7o) i s the 3-dimensional velocity dispersion of peaks identified on the filtering scale Rs 
(Eq. (|I5p. Hence, at the zeroth order, the evolved correlation of density peaks is simply obtained through a convolution 
of the initial correlation £ p k(i?s,^, r) with the diffusion kernel exp[— (l/3)fc 2 er^ pk (z)]. 



3. First order: linear growth due to coherent motions 



by 



At the first order, the Fourier transform of the covariance matrices E^ n '(q) receive a contribution 5E' Tn '(q) given 



1 



(B13) 



1 (e± • q) (e± • q) 2 



hpk(q) p 5s(<i) I i o 



o 1 



The velocity bias factor b vpk (q), Eq. (|16p . is the same for the two peak populations because it depends only on the 
filtering scale Rs- In addition, there is a non-zero mean velocity with line-of-sight components 



5~<?\ q ) = --^^-b vpk (q)b spk2 (q)P Ss (q) + ^J*— 



a-i q 

^ 0) (q) = — ^ 
<7_i q 

and transverse components 



b vpk (q)b spkl {q)Pg s (q) - 



2o-_i<t 2 
5i 
2ct_i(t 2 



3(r-q) 3 -(f -q)l b^q)^ P 5s (q) 



3 (f • q) 3 - (r • q)l b vpk (q)([ 0) P 5s (q) 



^ ±x) (q) 



^% (r • q) (e± • q) (e± • q) b vpk (q)(^P Ss (q) 



<J-\(J2 



^(q) = - 



— — g(r-q) (e± ■ q) (e± ■ q) ^rpk^CP^-^sC?) 



cr_icr 2 



(B14) 



(B15) 



Here, & S pki(<?) = + ^c^ 2 * s ^ ne linear bias for peaks of height v$ and curvature itj. Even if the peaks had all the 
same height, we would have to distinguish between the linear bias of the two populations because the peak curvatures 
would not necessarily be the same. We also remark that the mean transverse velocity components S^ 1 ^ vanish 
upon averaging over the orientations of the principal axis frames. With the aid of these results, we now expand the 
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right-hand side of Eq. (|B10l) up to first order and impose the peak constraint to arrive at 
^- y , d 6 Cirf 6 C2^ k (x' 1 )n pk (x^)P 2 (y 1 ,y 2 ;r')cxp (-^EJ + i^z) 



(2) 



n 



J- ld 6 C i n pk (xr)P 1 {y i 



(-ijt<JEJ + »jt<JH 



e -3 fc %1 1 



d 3 q 

(2tt) 3 



M<?)M<?) + — ' ' k ) M<z) + M?) fa vpk( 9 ) 



2 1.2 



c-x/ q 



(f • q) 2 (f' • k) 2 + J] (e' • q)(e' a • q)(e' a • k)« • k)l b 2 pk (<?) 



a=± 



P is (<?)e^ r (B16) 



upon an integration over the peak asymmetry parameters Vi, Wi and the peak curvatures Ui (see Appendix SjXj). In 
the second term of this equation, -Pi(yi) is a 1-point probability density and, in the third, f ■ e' ± = f ' ■ e' ± = where 
?' = r'/r'. It is important to bear in mind that the spectral moments a n are evaluated at initial redshift Zi, e.g. 
cr_i = cr_i(zj). Hence, the ratio a v /a^\ equals D{z) / D{zi). To perform the integral over r', we utilize the following 
relations 

/d 3 rf i f,e i ^- k > r = (27r) 3 ^ 3 )(q-k)fe4- 1 |rf 3 re ±1 e ±3 e !(q - k) ' r = (2^) 3 ^ 3 >(q - k)c ±! e ±3 , (B17) 

where e+ and e_ denote the unit vectors orthogonal to the wavevector k, i.e. k • e± = 0. As a consequence, all 
the terms involving e' ± • k or e' ± ■ k vanish. This cancellation reflects the fact that, on average, only streaming 
motions along the separation vector r' of a peak pair can affect the peak correlation £ p k- On employing the identity 
(f ■ q) 2 + X) a =±( e a ' <l)(®a ' 4) = 1 (which follows from Eq |A28[) . we thus obtain 



L. y"d 3 r'e- ikr ' Jd 6 ( 1 d 6 C2 n pk (^ 1 )n pk (x , 2 )P 2 (y 1 ,y 2 ;v') exp (-^EJ + *J f £ 



e -3 fe v bu(q)bi 2 (g) + 



t>n(g) + (512(9) bv P k(9) + 



g(g) 
D{zi) 



b 2 vpk (q)\Ps s (q). (B18) 



After some further manipulation, this yields the desired result Eq. (|63p where it is assumed v\ = v 2 = v 



4. Second order: Lagrangian and gravity mode-coupling 

The calculation of the 2nd order term in the Taylor expansion of Eq. (|B10|) , which is the lowest order contribution 
to the mode-coupling power, is long and fastidious. For clarity, we can decompose this second order mode-coupling 
power into three distinct pieces which we will compute successively : i) the second order contribution to the initial 
peak correlation £ p k(^i, v 2 , Rs, r ) convolved with the diffusion kernel, ii) (y\, v 2 , Rs, f) times the first order term 
in the expansion of exp(— ij^EJ + iJ^S) and hi) the second order term in the expansion of exp(— ^J^EJ + iJ'S). 

• The first piece reflects the fact that the initial, second order contribution is progressively smeared out by the peak 
motions. It is trivially 

J^e-^<^(u uV2 ,R s ,k)^ , (B19) 

(2) 

where P pk (vi, v 2 , Rg, k) is the Fourier transform of the second order correlation of initial density peaks, 
£ pk (vi,V2,Rs,r). 

• Taking the exponential damping factor out of the integral, the second part of the mode-coupling can be written 

as 

i=X,2 



^- /'d 6 Cin pk (xOP 1 (y i ) 
"pk J 



yjM- 1 BM- 1 y 2 ) --J f <5EJ + itfSZ 



(B20) 
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where the Fourier transform of y 2 M ! BM 1 yi is given by 



Jd 3 r (y+M-^M-Vi) e" iqr = j& spk i(# S pk 2 (g) - [3(r ■ q) 2 - l] (& spk i(?)d° } + W^Ki" 

+ [3(r. q ) 2 -l] 2 cl 0) d 0) + ^ 4 (f. q ) 2 

As we can see from Eq. ([B14p . there are also terms linear in and C^ 1 in the first order mean velocity 5E. On angle 
averaging over the variables which define the orientation of the principal axes, product of the form Ci°^Ci°^ reduce to 
(3/10) tt(Ci) as already shown in Appendix SjA] Hence, the Fourier transform of y|M^ 1 BM~ 1 yi and iJ^SE becomes 



(yjM^BM-Vi) (iJtJS)) = |(^7) (r-k)^^6 spkl (< 7l )6 spk2 (q 1 )[6 spkl ( g2 ) + 6 spk2 ( (Z2 )]fe vpk (q 2 
+ ^1 (£) I'Mr ■ k) [3(f • qi ) 2 - 1] [3(f • q 2 ) 3 - (f • q 2 )] 



5 S pkl 1 



(<Zi)tr(c 2 ) + ^(^^(C^jfovpkfe) >Ps a (qi)P6 3 fa) ■ ( B22 ) 



Upon substituting this result in Eq. (|B20p and integrating over the other variables, the second piece ii) can eventually 
be expressed as 



J(2n)*J (27r)3H -_ 1 ) q 2 2 



M ( f '. k )(£^) 

(7-1 J Q2 



(f ■ q 2 ) 2 (f ■ k) 2 + E « ■ ' ■ k)« • k) 



o=± 



fan(gi)bi 2 (gi)b 2 pk (g2) 



& S pki(gi)& sp ki(g2) bi 2 (gi) + 6 spk2 (gi)& spk2 (g 2 ) bn(gi) bvpkfe) 



1/(7. 



2o-s \(7_ 



» g 2 92 (f' • k) [3(f • qi ) 2 - 1] [3(f • q 2 ) 3 - (f ■ q 2 )] 



x I 9 Q lnG^ Q) (7i,7iz/ 2 )bii(gi) + 9 Q lnG^ Q) (7i,7ii/ 2 )bi2(gi) 



&v P kfe) ^ s (9i)^ s (g2)e t(qi+q2) - r ' , (B23) 



where we have replaced ti(Q) by (2/3)(3w 2 + w 2 ) before proceeding with the integration over the asymmetry param- 
eters (see Appendix [Af . Eq. (|A69|) can be employed to rewrite the average & P ki(9i)frpki(?2) in terms of the second 
order peak bias factors b 2 o, bn and bo 2 - 

• The third part of the mode-coupling power is the most difficult to compute. It can be written as 



e 3 K °v P k 



n 



" P k 



d 6 C i n pk (x;)P 1 (y < ) 



- f--jt<JEJ + irfdE J - ij f (5 2 EJ + irfd 2 E I . (B24) 



Let us first concentrate on the integral over the terms quadratic in <5E and fc. We must proceed carefully with the 
calculation of J^<5S because, once again, there are terms of the form Ci C; which do not vanish upon averaging 
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over the angular variables. A straightforward computation yields 



(Jt*3) 



(r-k) 



2 (r -qi) (f -q 2 ) 



45 / oy, 
2ct 2 i^ f7 _ 1 



,o=± 



tr 



(f • qx)(f ■ q 2 )(e + • k)(e_ • k) 
(Ci) +tr(c 2 2 )]6 vpk ( gi )6 vpk ( g2 ) \p 5s { qi )P 5s {q 2 ) 



&spkl(<7l) + frspk2(gi) &spkl(<72) + & S pk2(g2) &vpk(gi)&vpk(<?2) 

tr (Ci) + tr (cl) &vpk(gi)&v P k(g2) 

X! ( e « • ■ qi)( e -a ■ q2)(e_ a ' q2) 



(B25) 



Adding the contribution proportional to (J^EJ) 2 and (JtJEJ^Jtfc,) and integrating out the asymmetry parameters, 
we finally obtain 



'(2,)3i(2,)3H^J <&ll 



(r> ■ q 2 ) 2 (r' • k) 2 + J2 K ■ • q 2 )(e' a • k)(e' a ■ k) 



(f • qi ) 2 (f' • k) 2 + ^(e' • qi)(el • qj)(e' • k)(e' a ■ k) 

a=± 

K P M~bLM + 2(—) 3 K(i'-K 



a=± 



(f' • qi ) 2 (f • k) 2 + J2 « • qi)(e' a • qi)(e' a • k)(e' a • k) 



o=± 



<12 



biM+bM ^pk(9i)bvpk(g2)+( — ) (f'-k) 2 



jr' ■ qQ (f' 
<7i 



& sp ki(gi)&spki(g2) + &ii(gi)&i2(g2) + &ii(g2)&i2(gi)+& sp k2(gi)^spk2(g2)| & vpk(gi)&v P k(92)-^2 ( 

(J 2 V(T_1 



(f'-qi)(f'-q 2 ) 



(q iq2 (r' ■ k) 2 [3(r' • Ql ) 2 - l] [3(f • q 2 ) 2 - l] - I2(e' + ■ k)(e'_ • k) « • qi )(e' a • qi )(e'_ a • q 2 )(e'_ a • q 2 )] J 



x I 9 a lnG[, Q) (7i , 711/1) + a Q lnG^ ; (7i ,711*2) 



(«)/ 



bvpk(gi)b vp k(g2) ^ ( 5 s (9i)- p <5s(92)e 



i (qi+q2)-r' 



(B26) 



where, for the sake of completeness, we have included all the terms involving (e' ± • k) even though they will vanish 
when we carry out the integral over r'. Next, we consider the integral over the second order terms S 2 T, and <5 2 S. The 
Fourier transform of the second order contribution to the covariancc matrices are 



<5 2 S(°)( qi , q 2 ) = -L/-^(f ■ qi )(f • q 2 ) + ^ qi q 2 [3(f ■ qi ) 2 - l] [3(f ■ q 2 ) 2 - l] 

' 1 / v-i , 1 2 7i -1 
— (qiq2) + — qiq 2 q 1 1 



(B27) 



1 - fi\ 



x &v P k (qi )b vp k (q 2 ) Ps s (<?i ) Ps s (?2 ) I2 



OQ(T 2 



>(f • qi)(f • q 2 ) 



5 2 S( ±1 )(q 1 ,q 2 ) = 



3 15 

2 + — <7i<7 2 (r • qi)(f • q 2 ) 

°\ ° 2 



(e± • qi)(e± • qi)(e± • q 2 )(e± • q 2 ) 



x & vp k(gi)^vpk(g2)-P5 s (qi)Pss (92) h 



40 



whereas, for the line-of-sight components of the mean velocity, we arrive at 



cr-i 



1-7? 



-i / 1 



;9i9 2 



7i -i 2 7i 



0W2 



-91 q-2 



&0&2 



r g 2 (r-qi)(r-q 2 ) (B28) 



^2 9192 [3(r ' Qi) 2 - 1] [3(f • q 2 ) 2 - l] > (r • qi)6 sp ki(92)^ P k(9i)A s (9l)A s (92) 



^Hf Cm,*)- -H (1-7?) 

(T_i 1 



1/1-1 1 2 

+ T29192 



7i -1 2 7i 
"9i 9 2 - T^9i 



OW2 



&0&2 



r92(f • qi)(f • q 2 ) 



_ 4#2 qiq 2 [ 3 (r ' 4i) 2 - l] [3(r ■ q 2 ) 2 - l] \(r ■ qi)& sp k2(92)&v P k(9i)-P<5 s (gi)-P«5 s (92) 



We have ignored all terms linear in Q ^ because these will cancel out when we integrate over the angular variables. 

Furthermore, the second order contribution to the transverse velocity component, 5 2 Hj , can be ignored because 
it vanishes upon averaging over the orientations of the principal axis frames. On integrating over the asymmetry 
parameters and the peak curvature, the third piece iii) can eventually be cast into the form 



e -^ 2 < k f d ^ f d ( b 



(2tt) 3 J (2tt) 3 \a 
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[l-7{) (9192)- 1 - 
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2 ^2 



(f • qi)(f' • q 2 )(f' • k) 2 b V p k (gi)b vpk (g 2 ) 



2 - — 9i9a(f' • qi)(f ■ q 2 ) V (e' a ■ q.i)(e' a • qi)(e' t • q 2 )« • q 2 )(e^ • k)(e' t • k) b vpk ( g i)b vp k(92) 



0~» 
CT-1 
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A 52 ( f ' . q l)( r' . q 2 ) _ JL^ [ 3 (f . qi ) 2 - l] [3(f' • q 2 ) 2 - l] - (l - TiT* 



-if J_ , (gi92) 2 _ 7? 
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7 2 + 9l) 



(f • qi)(f ■ k) 2 b n (? 2 ) + b I2 (g 2 ) b vpk (9i) n s (9i)n s (9 2 )e i(qi+q2) ' r ' (B29) 



To derive the mode-coupling power Pmc(^), we must now add the contributions Eqs (|B23() . (|B26[) and (|B29[) and 
perform the integration over r'. At this point, it is convenient to express the results in terms of quantities at the 
collapse redshift zq rather than the initial redshift Zi 1. This change of fiducial redshift is readily achieved by 
making the replacement zi — > zq. 

Again, all the terms involving the multiplicative factors (e^_ ■ k) or (e^_ • k) cancel out and, as for the correlation of 
initial density peaks, the mode-coupling power can be drastically simplified upon substituting the expression of the 
second order peak bias 611(91,92,20), Eq. (|A69[) . Adding Eqs (|B23|) . (|B26|) and (|B29[) . the mode-coupling power can 
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be written 
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(B30) 



where, unless otherwise specified, all quantities are evaluated at redshift zq. In analogy with standard PT, we have 
defined the kernels T n as 



1 (k • qi) (k ■ q„) - - 
>n(qi, • • ■ q«) = — ; ; 6 vpk (9l) ■ ' ' 6 vpk (9„) 



n\ qi q n 
We have also introduced the symmetric function bj^ i (qi,q 2 ,z), 

1 ( D(zq) \ 



(B31) 



6m(9i,92,^) = -F 2 (qi,q 2 ) 



2 V D{z) J 



T\ (qi ) 6k (g 2 , z ) 4- T\ (q 2 ) 6ii (91 , z ) 



\ (^y) 2 6rii (91,92,*)), (B32) 



which represents the evolved (Eulerian), second order bias of initial density maxima in the Zel'dovich approximation. 
The action of 6^(91, g 2 , z) on fields and correlation functions is identical to that of 611(91, 92, zq) (see ^III A|) . With the 

aid of these auxiliary functions, we can rearrange the Fourier transform of (l^^Q^fanibm^g ' times the exponential 
damping with the first nine terms in the curly bracket of Eq. (|B30[) into the compact expression 

7A3 (73^) e-3 fe2 <*w y> qi y^^E^ ^ jZ) ^e^ _ (B33) 

Observing that —iq(r ■ q) and — (l/2)g 2 [3(f • q) 2 — 1] are the Fourier transform of ^j 1 ^ and £3 , the four last terms 
in the curly bracket of Eq. (|B30p can be combined in a similar way with some of the terms present in the initial peak 
correlation ^ p k(i?s,^, r) and eventually arrive at Eq. (|72j) . Upon changing to the variables x = qi/k and ji = k • qi, 
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the mode-coupling power spectrum can be explicitly written as 
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where the second order Eulerian peak bias is 
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2\x(l + x 2 -2x^i)\ L 'a 2 

x bio + boifc 2 (l + x 2 - 2xpi) 

'D( ZQ y 



1 - -§fc 2 (l + X 2 - 2^) 

D(z )\ (1-Xfi) 



D{Z )\ V( x _<± k 2 x 2 



D(z) ) x 



D(z) J (1 + x 2 - 2x/jl) 



1 - -§/s 2 (1 + x 2 - 2xfi) 



D(z) 



&20 + bnfc 2 (1 + 2x 2 - 2xn) + bo 2 k A x 2 (l + x 2 - 2xp) 



[bio + boik 2 x 2 ) 
(B35) 



To calculate the mode-coupling in configuration space, we simply Fourier transform Pmc(+> Rs> k, z ). 
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